WGCNA GO Analysis
Kristen Brown, modified ZD
1/29/2024
WGCNA GO Analysis
sessionInfo() #provides list of loaded packages and version of R. ## R version 4.3.2 (2023-10-31)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.0
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## loaded via a namespace (and not attached):
## [1] digest_0.6.34 R6_2.5.1 fastmap_1.1.1 xfun_0.41
## [5] cachem_1.0.8 knitr_1.45 htmltools_0.5.7 rmarkdown_2.25
## [9] lifecycle_1.0.4 cli_3.6.2 sass_0.4.8 jquerylib_0.1.4
## [13] compiler_4.3.2 rstudioapi_0.15.0 tools_4.3.2 evaluate_0.23
## [17] bslib_0.6.1 yaml_2.3.8 rlang_1.1.3 jsonlite_1.8.8First, load the necessary packages.
library(goseq)## Loading required package: BiasedUrn## Loading required package: geneLenDataBase## library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(forcats)
library(ggplot2)
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(tidyr)
library(grDevices)
library(reshape2)##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithslibrary(Rmisc)## Loading required package: lattice## Loading required package: plyr## ------------------------------------------------------------------------------## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)## ------------------------------------------------------------------------------##
## Attaching package: 'plyr'## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarizelibrary(ggpubr)##
## Attaching package: 'ggpubr'## The following object is masked from 'package:plyr':
##
## mutatelibrary(tibble)
library(gridExtra)
library(tidyr)
library(zoo)##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(ComplexHeatmap)## Loading required package: grid## ========================================
## ComplexHeatmap version 2.16.0
## Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
## Github page: https://github.com/jokergoo/ComplexHeatmap
## Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
##
## If you use it in published research, please cite either one:
## - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
## - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
## genomic data. Bioinformatics 2016.
##
##
## The new InteractiveComplexHeatmap package can directly export static
## complex heatmaps into an interactive Shiny app with zero effort. Have a try!
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(ComplexHeatmap))
## ========================================library(circlize)## ========================================
## circlize version 0.4.15
## CRAN page: https://cran.r-project.org/package=circlize
## Github page: https://github.com/jokergoo/circlize
## Documentation: https://jokergoo.github.io/circlize_book/book/
##
## If you use it in published research, please cite:
## Gu, Z. circlize implements and enhances circular visualization
## in R. Bioinformatics 2014.
##
## This message can be suppressed by:
## suppressPackageStartupMessages(library(circlize))
## ========================================library(GSEABase)## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following object is masked from 'package:gridExtra':
##
## combine## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.## Loading required package: annotate## Loading required package: AnnotationDbi## Loading required package: stats4## Loading required package: IRanges## Loading required package: S4Vectors##
## Attaching package: 'S4Vectors'## The following object is masked from 'package:plyr':
##
## rename## The following object is masked from 'package:tidyr':
##
## expand## The following objects are masked from 'package:dplyr':
##
## first, rename## The following object is masked from 'package:geneLenDataBase':
##
## unfactor## The following object is masked from 'package:utils':
##
## findMatches## The following objects are masked from 'package:base':
##
## expand.grid, I, unname##
## Attaching package: 'IRanges'## The following object is masked from 'package:plyr':
##
## desc## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice##
## Attaching package: 'AnnotationDbi'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: XML## Loading required package: graph##
## Attaching package: 'graph'## The following object is masked from 'package:XML':
##
## addNode## The following object is masked from 'package:circlize':
##
## degree## The following object is masked from 'package:plyr':
##
## joinlibrary(data.table)##
## Attaching package: 'data.table'## The following object is masked from 'package:IRanges':
##
## shift## The following objects are masked from 'package:S4Vectors':
##
## first, second## The following objects are masked from 'package:reshape2':
##
## dcast, melt## The following objects are masked from 'package:dplyr':
##
## between, first, lastlibrary(stringr)##
## Attaching package: 'stringr'## The following object is masked from 'package:graph':
##
## boundarylibrary(GenomicRanges)## Loading required package: GenomeInfoDblibrary(rtracklayer)
library(rrvgo)Load in data
library(rtracklayer)
gff<-rtracklayer::import("../../data/Pocillopora_acuta_HIv2.genes_fixed.gff3")
gff<-as.data.frame(gff)
dim(gff) # 478988 9## [1] 478988 13names(gff) ## [1] "seqnames" "start" "end" "width"
## [5] "strand" "source" "type" "score"
## [9] "phase" "ID" "transcript_id" "gene_id"
## [13] "Parent"transcripts <- subset(gff, type == "transcript")
transcripts_gr <- makeGRangesFromDataFrame(transcripts, keep.extra.columns=TRUE) #extract length information
transcript_lengths <- width(transcripts_gr) #isolate length of each gene
seqnames<-transcripts_gr$ID #extract list of gene id
lengths<-cbind(seqnames, transcript_lengths)
lengths<-as.data.frame(lengths) #convert to data frame
dim(transcripts) #33730 13## [1] 33730 13kegg <- read.delim("../../data/Pocillopora_acuta_HIv2.genes.KEGG_results.txt",header = FALSE)
kegg <- as.data.frame(kegg)
colnames(kegg)[1] <- "gene_id"
colnames(kegg)[2] <- "KEGG_new"
head(kegg)## gene_id KEGG_new
## 1 Pocillopora_acuta_HIv2___RNAseq.g24143.t1a
## 2 Pocillopora_acuta_HIv2___RNAseq.g24143.t1b K22584
## 3 Pocillopora_acuta_HIv2___RNAseq.g22918.t1
## 4 Pocillopora_acuta_HIv2___RNAseq.g18333.t1 K03386
## 5 Pocillopora_acuta_HIv2___RNAseq.g7985.t1
## 6 Pocillopora_acuta_HIv2___RNAseq.g13511.t1eggnog<-read.delim("../../data/Pocillopora_acuta_HIv2.genes.EggNog_results.txt")#this file contains all of the go terms, descriptions, kegg, etc
eggnog<- plyr::rename(eggnog, c("X.query"="gene_id"))
head(eggnog,2)## gene_id seed_ortholog evalue score
## 1 Pocillopora_acuta_HIv2___RNAseq.g24143.t1a 45351.EDO48725 2.16e-120 364
## 2 Pocillopora_acuta_HIv2___RNAseq.g18333.t1 45351.EDO38694 3.18e-123 355
## eggNOG_OGs
## 1 COG0620@1|root,KOG2263@2759|Eukaryota,38GZS@33154|Opisthokonta,3BNKS@33208|Metazoa
## 2 COG0450@1|root,KOG0852@2759|Eukaryota,38B9P@33154|Opisthokonta,3BGS4@33208|Metazoa
## max_annot_lvl COG_category
## 1 33208|Metazoa E
## 2 33208|Metazoa O
## Description Preferred_name
## 1 Cobalamin-independent synthase, Catalytic domain -
## 2 negative regulation of male germ cell proliferation PRDX4
## GOs
## 1 -
## 2 GO:0000003,GO:0001775,GO:0002252,GO:0002263,GO:0002274,GO:0002275,GO:0002283,GO:0002366,GO:0002376,GO:0002443,GO:0002444,GO:0002446,GO:0003006,GO:0003674,GO:0003824,GO:0004601,GO:0005488,GO:0005515,GO:0005575,GO:0005576,GO:0005615,GO:0005622,GO:0005623,GO:0005737,GO:0005783,GO:0005790,GO:0005829,GO:0006082,GO:0006457,GO:0006464,GO:0006468,GO:0006520,GO:0006575,GO:0006793,GO:0006796,GO:0006807,GO:0006810,GO:0006887,GO:0006915,GO:0006950,GO:0006952,GO:0006955,GO:0006979,GO:0007154,GO:0007165,GO:0007249,GO:0007252,GO:0007275,GO:0007276,GO:0007283,GO:0007548,GO:0008150,GO:0008152,GO:0008219,GO:0008285,GO:0008379,GO:0008406,GO:0008584,GO:0009056,GO:0009266,GO:0009409,GO:0009605,GO:0009607,GO:0009617,GO:0009628,GO:0009636,GO:0009893,GO:0009966,GO:0009967,GO:0009987,GO:0010467,GO:0010604,GO:0010646,GO:0010647,GO:0010941,GO:0010942,GO:0010950,GO:0010952,GO:0012501,GO:0012505,GO:0016043,GO:0016192,GO:0016209,GO:0016310,GO:0016491,GO:0016684,GO:0016999,GO:0017001,GO:0017144,GO:0019222,GO:0019471,GO:0019538,GO:0019725,GO:0019752,GO:0019953,GO:0022414,GO:0022417,GO:0023051,GO:0023052,GO:0023056,GO:0030141,GO:0030162,GO:0030198,GO:0031323,GO:0031325,GO:0031410,GO:0031974,GO:0031982,GO:0031983,GO:0032268,GO:0032270,GO:0032501,GO:0032502,GO:0032504,GO:0032940,GO:0033554,GO:0034774,GO:0035556,GO:0036211,GO:0036230,GO:0042119,GO:0042127,GO:0042221,GO:0042592,GO:0042737,GO:0042742,GO:0042743,GO:0042744,GO:0042802,GO:0042803,GO:0042981,GO:0043062,GO:0043065,GO:0043067,GO:0043068,GO:0043085,GO:0043170,GO:0043207,GO:0043226,GO:0043227,GO:0043229,GO:0043231,GO:0043233,GO:0043280,GO:0043281,GO:0043299,GO:0043312,GO:0043412,GO:0043436,GO:0043900,GO:0043901,GO:0044093,GO:0044237,GO:0044238,GO:0044248,GO:0044260,GO:0044267,GO:0044281,GO:0044421,GO:0044422,GO:0044424,GO:0044433,GO:0044444,GO:0044446,GO:0044464,GO:0044703,GO:0045055,GO:0045137,GO:0045321,GO:0045454,GO:0045862,GO:0046425,GO:0046427,GO:0046546,GO:0046661,GO:0046903,GO:0046983,GO:0048232,GO:0048513,GO:0048518,GO:0048519,GO:0048522,GO:0048523,GO:0048583,GO:0048584,GO:0048608,GO:0048609,GO:0048731,GO:0048856,GO:0050789,GO:0050790,GO:0050794,GO:0050896,GO:0051171,GO:0051173,GO:0051179,GO:0051186,GO:0051187,GO:0051234,GO:0051239,GO:0051241,GO:0051246,GO:0051247,GO:0051336,GO:0051345,GO:0051604,GO:0051704,GO:0051707,GO:0051716,GO:0051920,GO:0052547,GO:0052548,GO:0055114,GO:0060205,GO:0060255,GO:0061458,GO:0065007,GO:0065008,GO:0065009,GO:0070013,GO:0070417,GO:0070887,GO:0071704,GO:0071840,GO:0072593,GO:0080090,GO:0097190,GO:0097237,GO:0097708,GO:0098542,GO:0098754,GO:0098869,GO:0099503,GO:0101002,GO:1901564,GO:1901605,GO:1902531,GO:1902533,GO:1904813,GO:1904892,GO:1904894,GO:1905936,GO:1905937,GO:1990748,GO:2000116,GO:2000241,GO:2000242,GO:2000254,GO:2000255,GO:2001056,GO:2001233,GO:2001235,GO:2001267,GO:2001269
## EC KEGG_ko
## 1 2.1.1.14 ko:K00549
## 2 1.11.1.15 ko:K03386
## KEGG_Pathway
## 1 ko00270,ko00450,ko01100,ko01110,ko01230,map00270,map00450,map01100,map01110,map01230
## 2 ko04214,map04214
## KEGG_Module KEGG_Reaction KEGG_rclass
## 1 M00017 R04405,R09365 RC00035,RC00113,RC01241
## 2 - - -
## BRITE KEGG_TC CAZy BiGG_Reaction
## 1 ko00000,ko00001,ko00002,ko01000 - - -
## 2 ko00000,ko00001,ko01000,ko04147 - - -
## PFAMs
## 1 Meth_synt_2
## 2 1-cysPrx_C,AhpC-TSAgogene <- merge(transcripts, eggnog, by=c("gene_id"), all=T)
gogene <- merge(gogene, kegg, by=c("gene_id"), all=T)
head(gogene,2)## gene_id seqnames
## 1 Pocillopora_acuta_HIv2___RNAseq.10002_t Pocillopora_acuta_HIv2___Sc0000013
## 2 Pocillopora_acuta_HIv2___RNAseq.10010_t Pocillopora_acuta_HIv2___Sc0000013
## start end width strand source type score.x phase
## 1 4542087 4551503 9417 + AUGUSTUS transcript NA NA
## 2 4639103 4647350 8248 + AUGUSTUS transcript NA NA
## ID
## 1 Pocillopora_acuta_HIv2___RNAseq.10002_t
## 2 Pocillopora_acuta_HIv2___RNAseq.10010_t
## transcript_id Parent seed_ortholog evalue
## 1 Pocillopora_acuta_HIv2___RNAseq.10002_t 45351.EDO27354 2.41e-93
## 2 Pocillopora_acuta_HIv2___RNAseq.10010_t 6087.XP_002166004.2 1.28e-38
## score.y
## 1 317
## 2 164
## eggNOG_OGs
## 1 COG0666@1|root,KOG0510@2759|Eukaryota,38G7Q@33154|Opisthokonta,3BCDU@33208|Metazoa
## 2 COG0666@1|root,KOG4177@2759|Eukaryota
## max_annot_lvl COG_category Description
## 1 33208|Metazoa DZ osmolarity-sensing cation channel activity
## 2 2759|Eukaryota I spectrin binding
## Preferred_name
## 1 TRPA1
## 2 -
## GOs
## 1 GO:0000302,GO:0001580,GO:0002791,GO:0002793,GO:0003008,GO:0003012,GO:0003674,GO:0004888,GO:0005034,GO:0005215,GO:0005216,GO:0005217,GO:0005244,GO:0005245,GO:0005261,GO:0005262,GO:0005488,GO:0005515,GO:0005575,GO:0005623,GO:0005886,GO:0005887,GO:0006810,GO:0006811,GO:0006812,GO:0006816,GO:0006873,GO:0006874,GO:0006875,GO:0006936,GO:0006939,GO:0006950,GO:0006979,GO:0007154,GO:0007165,GO:0007166,GO:0007204,GO:0007600,GO:0007602,GO:0007606,GO:0007610,GO:0007638,GO:0008150,GO:0008324,GO:0009266,GO:0009314,GO:0009408,GO:0009409,GO:0009410,GO:0009416,GO:0009453,GO:0009581,GO:0009582,GO:0009583,GO:0009593,GO:0009605,GO:0009612,GO:0009628,GO:0009636,GO:0009719,GO:0009966,GO:0009967,GO:0009987,GO:0010033,GO:0010035,GO:0010037,GO:0010243,GO:0010378,GO:0010646,GO:0010647,GO:0010817,GO:0014070,GO:0014074,GO:0014832,GO:0014848,GO:0015075,GO:0015085,GO:0015267,GO:0015276,GO:0015278,GO:0015318,GO:0016020,GO:0016021,GO:0016043,GO:0016048,GO:0016324,GO:0019233,GO:0019722,GO:0019725,GO:0019932,GO:0022607,GO:0022803,GO:0022832,GO:0022834,GO:0022836,GO:0022838,GO:0022839,GO:0022843,GO:0022857,GO:0022890,GO:0023041,GO:0023051,GO:0023052,GO:0023056,GO:0030001,GO:0030003,GO:0030424,GO:0031000,GO:0031224,GO:0031226,GO:0031644,GO:0031646,GO:0032024,GO:0032421,GO:0032501,GO:0032879,GO:0032880,GO:0032991,GO:0033554,GO:0033555,GO:0034220,GO:0034605,GO:0034702,GO:0034703,GO:0035556,GO:0035690,GO:0035774,GO:0036270,GO:0038023,GO:0040011,GO:0040040,GO:0042221,GO:0042330,GO:0042331,GO:0042391,GO:0042493,GO:0042542,GO:0042592,GO:0042752,GO:0042802,GO:0042995,GO:0043005,GO:0043052,GO:0043269,GO:0043270,GO:0043279,GO:0043933,GO:0044057,GO:0044070,GO:0044085,GO:0044425,GO:0044459,GO:0044464,GO:0045177,GO:0046677,GO:0046873,GO:0046883,GO:0046887,GO:0046957,GO:0048265,GO:0048518,GO:0048519,GO:0048522,GO:0048523,GO:0048583,GO:0048584,GO:0048878,GO:0050708,GO:0050714,GO:0050789,GO:0050794,GO:0050796,GO:0050801,GO:0050848,GO:0050850,GO:0050877,GO:0050896,GO:0050906,GO:0050907,GO:0050909,GO:0050912,GO:0050913,GO:0050951,GO:0050954,GO:0050955,GO:0050960,GO:0050961,GO:0050965,GO:0050966,GO:0050968,GO:0050974,GO:0050982,GO:0051046,GO:0051047,GO:0051049,GO:0051050,GO:0051179,GO:0051209,GO:0051222,GO:0051223,GO:0051234,GO:0051239,GO:0051240,GO:0051259,GO:0051260,GO:0051262,GO:0051282,GO:0051283,GO:0051289,GO:0051480,GO:0051606,GO:0051641,GO:0051649,GO:0051716,GO:0051930,GO:0051931,GO:0051969,GO:0052129,GO:0055065,GO:0055074,GO:0055080,GO:0055082,GO:0055085,GO:0060089,GO:0060341,GO:0060401,GO:0060402,GO:0061178,GO:0065003,GO:0065007,GO:0065008,GO:0070201,GO:0070417,GO:0070588,GO:0070838,GO:0070887,GO:0071241,GO:0071244,GO:0071310,GO:0071312,GO:0071313,GO:0071407,GO:0071415,GO:0071417,GO:0071466,GO:0071495,GO:0071840,GO:0071944,GO:0072347,GO:0072503,GO:0072507,GO:0072511,GO:0090087,GO:0090276,GO:0090277,GO:0097458,GO:0097553,GO:0097603,GO:0097604,GO:0098590,GO:0098655,GO:0098660,GO:0098662,GO:0098771,GO:0098796,GO:0098862,GO:0098900,GO:0098908,GO:0099094,GO:0099604,GO:0120025,GO:1901698,GO:1901699,GO:1901700,GO:1901701,GO:1902495,GO:1902531,GO:1902533,GO:1903522,GO:1903530,GO:1903532,GO:1903793,GO:1904058,GO:1904951,GO:1990351,GO:1990760
## 2 -
## EC KEGG_ko KEGG_Pathway KEGG_Module KEGG_Reaction KEGG_rclass
## 1 - ko:K04984 ko04750,map04750 - - -
## 2 - - - - - -
## BRITE KEGG_TC CAZy
## 1 ko00000,ko00001,ko04040 1.A.4.6.1,1.A.4.6.2,1.A.4.6.3,1.A.4.6.5 -
## 2 - - -
## BiGG_Reaction PFAMs KEGG_new
## 1 - Ank,Ank_2,Ank_3,Ank_4,Ion_trans K04984
## 2 - Ank,Ank_2,Ank_4,Ank_5,Ion_trans K04984geneInfo <- read.csv("../../output/WGCNA/WGCNA_ModuleMembership.csv") #this file was generated from the WGCNA analyses and contains the modules of interest
geneInfo<- plyr::rename(geneInfo, c("X"="gene_id"))
dim(geneInfo) # there are 9012 genes in our gene info file## [1] 9012 40geneInfo <- merge(gogene, geneInfo, by=c("gene_id"), all=T)Format GO terms to remove dashes and quotes and separate by semicolons (replace , with ;) in GOs column
geneInfo$GOs <- gsub(",", ";", geneInfo$GOs)
geneInfo$GOs <- gsub('"', "", geneInfo$GOs)
geneInfo$GOs <- gsub("-", NA, geneInfo$GOs)
geneInfo$KEGG_new[geneInfo$KEGG_new == ""] <- NAunique(geneInfo$moduleColor)## [1] "green" NA "blue" "salmon" "turquoise"
## [6] "yellow" "black" "red" "magenta" "lightcyan"
## [11] "purple" "brown" "pink" "midnightblue" "tan"
## [16] "cyan"geneInfo$Length<-lengths$transcript_lengths[match(geneInfo$gene_id, lengths$seqnames)]Run the GOSeq function by module color to test for GO term
enrichment. Due to high number of enriched GO terms, I am using an
adjusted p-value threshold of <0.001 and using rrvgo
package to reduce redundancy in list of GO terms.
Calcification
### Generate vector with names of all genes
ALL.vector <- c(geneInfo$gene_id)### Generate length vector for all genes
LENGTH.vector <- as.integer(geneInfo$Length)calc_up_mods <- c("brown", "red", "black", "pink", "salmon", "blue")
nrow(geneInfo %>% dplyr::filter(moduleColor=="brown")) #942## [1] 942nrow(geneInfo %>% dplyr::filter(moduleColor=="red")) #425## [1] 425nrow(geneInfo %>% filter(moduleColor=="black")) #396## [1] 396nrow(geneInfo %>% filter(moduleColor=="pink")) #220## [1] 220nrow(geneInfo %>% filter(moduleColor=="salmon")) #154## [1] 154nrow(geneInfo %>% filter(moduleColor=="blue")) #1989## [1] 1989sum(nrow(geneInfo %>% dplyr::filter(moduleColor=="brown")), nrow(geneInfo %>% dplyr::filter(moduleColor=="red")), nrow(geneInfo %>% filter(moduleColor=="black")), nrow(geneInfo %>% filter(moduleColor=="pink")), nrow(geneInfo %>% filter(moduleColor=="salmon")), nrow(geneInfo %>% filter(moduleColor=="blue")))## [1] 4126# 4126
calc_down_mods <- c("turquoise","magenta","lightcyan")
nrow(geneInfo %>% dplyr::filter(moduleColor=="turquoise")) #2558## [1] 2558nrow(geneInfo %>% dplyr::filter(moduleColor=="magenta")) #219## [1] 219nrow(geneInfo %>% filter(moduleColor=="lightcyan")) #65## [1] 65sum(nrow(geneInfo %>% dplyr::filter(moduleColor=="turquoise")), nrow(geneInfo %>% dplyr::filter(moduleColor=="magenta")), nrow(geneInfo %>% filter(moduleColor=="lightcyan")))## [1] 2842# 2842
other_mods <- c("green","yellow", "purple", "midnightblue","cyan","tan")
sum(nrow(geneInfo %>% dplyr::filter(moduleColor=="green")), nrow(geneInfo %>% dplyr::filter(moduleColor=="yellow")), nrow(geneInfo %>% filter(moduleColor=="purple")), nrow(geneInfo %>% filter(moduleColor=="midnightblue")), nrow(geneInfo %>% filter(moduleColor=="cyan")),nrow(geneInfo %>% filter(moduleColor=="tan")))## [1] 2044# 2044
# 4126 + 2842 + 2044 = 9012, which represents all of our genes4126 genes are in the 6 modules significantly upregulated by calcification.
### Generate vector with names in just the module we are analyzing
# ID.vector <- geneInfo %>%
# filter(moduleColor==c("brown", "red", "black", "pink", "salmon", "green")) %>%
# #get_rows(.data[[module]]))%>%
# pull(gene_id)
ID.vector <- geneInfo %>%
filter(moduleColor %in% c("brown", "red", "black", "pink", "salmon", "blue")) %>%
pull(gene_id)
length(ID.vector) #4126## [1] 4126##Get a list of GO Terms for each module
GO.terms <- geneInfo %>%
filter(moduleColor %in% c("brown", "red", "black", "pink", "salmon", "blue")) %>%
#filter(get_rows(.data[[module]]))%>%
dplyr::select(GOs,gene_id) %>% rename(GOs = "GO.terms")
dim(GO.terms) #4126 2## [1] 4126 2##Format to have one goterm per row with gene ID repeated
split <- strsplit(as.character(GO.terms$GO.terms), ";")
split2 <- data.frame(v1 = rep.int(GO.terms$gene, sapply(split, length)), v2 = unlist(split))
colnames(split2) <- c("gene", "GO.terms")
GO.terms<-split2##Construct list of genes with 1 for genes in module and 0 for genes not in the module
gene.vector=as.integer(ALL.vector %in% ID.vector)
names(gene.vector)<-ALL.vector#set names
#weight gene vector by bias for length of gene
pwf<-nullp(gene.vector, ID.vector, bias.data=LENGTH.vector) ## Warning in pcls(G): initial point very close to some inequality constraints
#run goseq using Wallenius method for all categories of GO terms
GO.wall<-goseq(pwf, ID.vector, gene2cat=GO.terms, test.cats=c("GO:BP", "GO:MF", "GO:CC"), method="Wallenius", use_genes_without_cat=TRUE)## Using manually entered categories.## Calculating the p-values...## 'select()' returned 1:1 mapping between keys and columnsGO <- GO.wall[order(GO.wall$over_represented_pvalue),]
colnames(GO)[1] <- "GOterm"#adjust p-values
GO$bh_adjust <- p.adjust(GO$over_represented_pvalue, method="BH") #add adjusted p-values#Filtering for p < 0.01
GO <- GO %>%
dplyr::filter(bh_adjust<0.00001) %>%
dplyr::arrange(., ontology, bh_adjust)
#Write file of results
write.csv(GO, file = "../../output/WGCNA/GO_analysis/goseq_pattern_calcification.csv")#add vector for terms of interest to reduce number of GO terms - NOT using this to look at individual modules for exploratory purposes
keywords<-c("metabolism", "carbon","bicarbonate", "apoptosis", "death", "symbiosis", "regulation of cell communication", "trans membrane transport", "transmembrane", "organic substance transport", "inorganic substance transport","response to stress", "antioxidant", "calcification","biomineralization", "heat","transporters","proton transport","ion transport","acid-base", "homeostasis")go_results <- read.csv("../../output/WGCNA/GO_analysis/goseq_pattern_calcification.csv") go_results<-go_results%>%
filter(ontology=="BP")%>%
filter(bh_adjust != "NA") %>%
filter(numInCat>10)%>%
arrange(., bh_adjust)
dim(go_results)## [1] 2774 9head(go_results)## X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 GO:0000003 0 1 465
## 2 2 GO:0006139 0 1 608
## 3 3 GO:0006355 0 1 462
## 4 4 GO:0006725 0 1 650
## 5 5 GO:0006807 0 1 1215
## 6 6 GO:0006810 0 1 664
## numInCat term ontology bh_adjust
## 1 465 reproduction BP 0
## 2 608 nucleobase-containing compound metabolic process BP 0
## 3 462 regulation of DNA-templated transcription BP 0
## 4 650 cellular aromatic compound metabolic process BP 0
## 5 1215 nitrogen compound metabolic process BP 0
## 6 664 transport BP 0library(rrvgo)
#Reduce/collapse GO term set with the rrvgo package
simMatrix <- calculateSimMatrix(go_results$GOterm,
orgdb="org.Ce.eg.db", #c. elegans database
ont="BP",
method="Rel")## ## preparing gene to GO mapping data...## preparing IC data... #calculate similarity
scores <- setNames(-log(go_results$bh_adjust), go_results$GOterm)
reducedTerms <- reduceSimMatrix(simMatrix,
scores,
threshold=0.7,
orgdb="org.Ce.eg.db")
dim(reducedTerms)## [1] 2213 10#keep only the goterms from the reduced list
go_results<-go_results%>%
filter(GOterm %in% reducedTerms$go)
#add in parent terms to list of go terms
go_results$ParentTerm<-reducedTerms$parentTerm[match(go_results$GOterm, reducedTerms$go)]write.csv(go_results, "../../output/WGCNA/GO_analysis/goseq_pattern_calcification_filtered.csv")#plot significantly enriched GO terms by Slim Category faceted by slim term
GO.plot_calcification <- ggplot(go_results, aes(x = ontology, y = term)) +
geom_point(aes(size=bh_adjust)) +
scale_size(name="Over rep. p-value", trans="reverse", range=c(1,3))+
facet_grid(ParentTerm ~ ., scales = "free", labeller = label_wrap_gen(width = 5, multi_line = TRUE))+
theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
strip.text.y = element_text(angle=0, size = 10),
strip.text.x = element_text(size = 20),
axis.text = element_text(size = 8),
axis.title.x = element_blank(),
axis.title.y = element_blank())
GO.plot_calcification
ggsave(filename="../../output/WGCNA/GO_analysis/GO.plot_calcification.png", plot=GO.plot_calcification, dpi=300, height=100, units="in", limitsize=FALSE)## Saving 7 x 100 in imageCount number of GOterms by ParentTerm for the upregulation of calcification
library(dplyr)
result <- go_results %>%
dplyr::group_by(ParentTerm) %>%
dplyr::summarize(Number.of.terms = n_distinct(term))%>%
mutate(Calcification.direction = "Up")
print(result)## # A tibble: 149 × 3
## ParentTerm Number.of.terms Calcification.direct…¹
## <chr> <int> <chr>
## 1 DNA metabolic process 26 Up
## 2 DNA-templated transcription initiation 6 Up
## 3 RNA processing 12 Up
## 4 actin filament-based process 12 Up
## 5 aging 2 Up
## 6 amide metabolic process 12 Up
## 7 ammonium ion metabolic process 1 Up
## 8 anatomical structure morphogenesis 27 Up
## 9 animal organ development 15 Up
## 10 behavior 10 Up
## # ℹ 139 more rows
## # ℹ abbreviated name: ¹​Calcification.directionCalcification - down
2842 genes are in the 4 modules significantly downregulated by calcification.
### Generate vector with names of all genes
ALL.vector <- c(geneInfo$gene_id)### Generate length vector for all genes
LENGTH.vector <- as.integer(geneInfo$Length)### Generate vector with names in just the module we are analyzing
ID.vector <- geneInfo%>%
filter(moduleColor %in% c("turquoise","magenta","lightcyan"))%>%
pull(gene_id)
length(ID.vector) #2842## [1] 2842##Get a list of GO Terms for each module
GO.terms <- geneInfo%>%
filter(moduleColor %in% c("turquoise","magenta","lightcyan"))%>%
dplyr::select(GOs,gene_id) %>% rename(GOs = "GO.terms")
dim(GO.terms) #2842 2## [1] 2842 2##Format to have one goterm per row with gene ID repeated
split <- strsplit(as.character(GO.terms$GO.terms), ";")
split2 <- data.frame(v1 = rep.int(GO.terms$gene, sapply(split, length)), v2 = unlist(split))
colnames(split2) <- c("gene", "GO.terms")
GO.terms<-split2##Construct list of genes with 1 for genes in module and 0 for genes not in the module
gene.vector=as.integer(ALL.vector %in% ID.vector)
names(gene.vector)<-ALL.vector#set names
#weight gene vector by bias for length of gene
pwf<-nullp(gene.vector, ID.vector, bias.data=LENGTH.vector) ## Warning in pcls(G): initial point very close to some inequality constraints
#run goseq using Wallenius method for all categories of GO terms
GO.wall<-goseq(pwf, ID.vector, gene2cat=GO.terms, test.cats=c("GO:BP", "GO:MF", "GO:CC"), method="Wallenius", use_genes_without_cat=TRUE)## Using manually entered categories.## Calculating the p-values...## 'select()' returned 1:1 mapping between keys and columnsGO <- GO.wall[order(GO.wall$over_represented_pvalue),]
colnames(GO)[1] <- "GOterm"#adjust p-values
GO$bh_adjust <- p.adjust(GO$over_represented_pvalue, method="BH") #add adjusted p-values#Filtering for p < 0.01
GO <- GO %>%
dplyr::filter(bh_adjust<0.00001) %>%
dplyr::arrange(., ontology, bh_adjust)
#Write file of results
write.csv(GO, file = "../../output/WGCNA/GO_analysis/goseq_pattern_calcification_down.csv")go_results <-read.csv("../../output/WGCNA/GO_analysis/goseq_pattern_calcification_down.csv")go_results<-go_results%>%
filter(ontology=="BP")%>%
filter(bh_adjust != "NA") %>%
filter(numInCat>10)%>%
arrange(., bh_adjust)
dim(go_results)## [1] 1892 9head(go_results)## X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 GO:0006139 0 1 485
## 2 2 GO:0006725 0 1 513
## 3 3 GO:0006807 0 1 910
## 4 4 GO:0006810 0 1 419
## 5 5 GO:0006950 0 1 364
## 6 6 GO:0006996 0 1 461
## numInCat term ontology bh_adjust
## 1 485 nucleobase-containing compound metabolic process BP 0
## 2 513 cellular aromatic compound metabolic process BP 0
## 3 910 nitrogen compound metabolic process BP 0
## 4 419 transport BP 0
## 5 364 response to stress BP 0
## 6 461 organelle organization BP 0library(rrvgo)
#Reduce/collapse GO term set with the rrvgo package
simMatrix <- calculateSimMatrix(go_results$GOterm,
orgdb="org.Ce.eg.db", #c. elegans database
ont="BP",
method="Rel")## preparing gene to GO mapping data...## preparing IC data... #calculate similarity
scores <- setNames(-log(go_results$bh_adjust), go_results$GOterm)
reducedTerms <- reduceSimMatrix(simMatrix,
scores,
threshold=0.7,
orgdb="org.Ce.eg.db")
dim(reducedTerms)## [1] 1638 10#keep only the goterms from the reduced list
go_results<-go_results%>%
filter(GOterm %in% reducedTerms$go)
#add in parent terms to list of go terms
go_results$ParentTerm<-reducedTerms$parentTerm[match(go_results$GOterm, reducedTerms$go)]write.csv(go_results, "../../output/WGCNA/GO_analysis/goseq_pattern_calcification_down_filtered.csv")#plot significantly enriched GO terms by Slim Category faceted by slim term
GO.plot_calcification_down <- ggplot(go_results, aes(x = ontology, y = term)) +
geom_point(aes(size=bh_adjust)) +
scale_size(name="Over rep. p-value", trans="reverse", range=c(1,3))+
facet_grid(ParentTerm ~ ., scales = "free", labeller = label_wrap_gen(width = 5, multi_line = TRUE))+
theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
strip.text.y = element_text(angle=0, size = 10),
strip.text.x = element_text(size = 20),
axis.text = element_text(size = 8),
axis.title.x = element_blank(),
axis.title.y = element_blank())
GO.plot_calcification_down
ggsave(filename="../../output/WGCNA/GO_analysis/GO.plot_calcification_down.png", plot=GO.plot_calcification_down, dpi=300, height=100, units="in", limitsize=FALSE)## Saving 7 x 100 in imageCount number of GOterms by ParentTerm for the downregulation of calcification
library(dplyr)
result_down <- go_results %>%
dplyr::group_by(ParentTerm) %>%
dplyr::summarize(Number.of.terms = n_distinct(term))%>%
mutate(Calcification.direction = "Down")
print(result_down)## # A tibble: 128 × 3
## ParentTerm Number.of.terms Calcification.direct…¹
## <chr> <int> <chr>
## 1 RNA processing 29 Down
## 2 actin filament-based process 11 Down
## 3 aging 2 Down
## 4 amide metabolic process 12 Down
## 5 ammonium ion metabolic process 1 Down
## 6 anatomical structure morphogenesis 29 Down
## 7 animal organ development 16 Down
## 8 behavior 8 Down
## 9 biological process involved in inters… 7 Down
## 10 biosynthetic process 24 Down
## # ℹ 118 more rows
## # ℹ abbreviated name: ¹​Calcification.directioncode for by module
library(rrvgo)
# Define the unique module colors
module_colors <- na.omit(unique(geneInfo$moduleColor))
# Generate vector with names of all genes
ALL.vector <- c(geneInfo$gene_id)
# Generate length vector for all genes
LENGTH.vector <- as.integer(geneInfo$Length)
# Loop over each unique module color
for (color in module_colors) {
# Filter geneInfo based on the current color
color_filtered <- geneInfo %>% filter(moduleColor == color)
# Generate vector with names in just the module we are analyzing
ID.vector <- color_filtered$gene_id
length(ID.vector)
# Get a list of GO Terms for each module
GO.terms <- color_filtered %>%
dplyr::select(GOs, gene_id) %>%
rename(GOs = "GO.terms")
dim(GO.terms)
## Format to have one GO term per row with gene ID repeated
split <- strsplit(as.character(GO.terms$GO.terms), ";")
split2 <- data.frame(v1 = rep.int(GO.terms$gene_id, sapply(split, length)), v2 = unlist(split))
colnames(split2) <- c("gene", "GO.terms")
GO.terms <- split2
## Construct list of genes with 1 for genes in module and 0 for genes not in the module
gene.vector <- as.integer(ALL.vector %in% ID.vector)
names(gene.vector) <- ALL.vector # set names
# Weight gene vector by bias for length of gene
pwf <- nullp(gene.vector, ID.vector, bias.data = LENGTH.vector)
# Run goseq using Wallenius method for all categories of GO terms
GO.wall <- goseq(pwf, ID.vector, gene2cat = GO.terms, test.cats = c("GO:BP", "GO:MF", "GO:CC"), method = "Wallenius", use_genes_without_cat = TRUE)
GO <- GO.wall[order(GO.wall$over_represented_pvalue),]
colnames(GO)[1] <- "GOterm"
# Adjust p-values
GO$bh_adjust <- p.adjust(GO$over_represented_pvalue, method = "BH")
# Filtering for p < 0.01
GO <- GO %>%
dplyr::filter(bh_adjust < 0.00001) %>%
dplyr::arrange(., ontology, bh_adjust)
# Write file of results
write.csv(GO, file = paste0("../../output/WGCNA/GO_analysis/goseq_pattern_", color, ".csv"))
go_results <- GO
go_results<-go_results%>%
filter(ontology=="BP")%>%
filter(bh_adjust != "NA") %>%
filter(numInCat>100)%>%
arrange(., bh_adjust)
#Reduce/collapse GO term set with the rrvgo package
simMatrix <- calculateSimMatrix(go_results$GOterm,
orgdb="org.Ce.eg.db", #c. elegans database
ont="BP",
method="Rel")
#calculate similarity
scores <- setNames(-log(go_results$bh_adjust), go_results$GOterm)
reducedTerms <- reduceSimMatrix(simMatrix,
scores,
threshold=0.7,
orgdb="org.Ce.eg.db")
#keep only the goterms from the reduced list
go_results <- go_results %>% filter(GOterm %in% reducedTerms$go)
#add in parent terms to list of go terms
go_results$ParentTerm<-reducedTerms$parentTerm[match(go_results$GOterm, reducedTerms$go)]
#plot significantly enriched GO terms by Slim Category faceted by slim term
GO.plot <- ggplot(go_results, aes(x = ontology, y = term)) +
geom_point(aes(size=bh_adjust)) +
scale_size(name="Over rep. p-value", trans="reverse", range=c(1,3))+
facet_grid(ParentTerm ~ ., scales = "free", labeller = label_wrap_gen(width = 5, multi_line = TRUE))+
theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
strip.text.y = element_text(angle=0, size = 10),
strip.text.x = element_text(size = 20),
axis.text = element_text(size = 8),
axis.title.x = element_blank(),
axis.title.y = element_blank())
GO.plot
length(colnames(go_results)[go_results$ParentTerm=="cation transport"])
length(colnames(go_results)[go_results$ParentTerm=="inorganic ion homeostasis"])
length(colnames(go_results)[go_results$ParentTerm=="regulation of cellular response to stress"])
}#Specific differentially expressed genes
wgcna_counts_filtered<-read.csv("../../output/Filtered_gene_count_matrix.csv", strip.white=T)
wgcna_counts_filtered<- plyr::rename(wgcna_counts_filtered, c("X"="Gene"))colnames(wgcna_counts_filtered)## [1] "Gene" "RF13B" "RF13D" "RF14B" "RF14C" "RF15B" "RF15D" "RF17B" "RF17D"
## [10] "RF18B" "RF18D" "RF19B" "RF19C" "RF20B" "RF20C" "RF22B" "RF22C" "RF23A"
## [19] "RF23C" "RF24B" "RF24D" "RF25A" "RF25C" "RS11B" "RS11D" "RS12A" "RS12C"
## [28] "RS13A" "RS13C" "RS14B" "RS14C" "RS15B" "RS15D" "RS1B" "RS1C" "RS2B"
## [37] "RS2C" "RS3B" "RS3D" "RS6A" "RS6D" "RS7B" "RS7C" "RS8B" "RS8C"
## [46] "RS9A" "RS9C"library(tidyr)
wgcna_counts_filtered_long <- pivot_longer(wgcna_counts_filtered, cols=2:47, names_to = "Colony", values_to = "Counts")
wgcna_counts_filtered_long$Colony <- as.factor(wgcna_counts_filtered_long$Colony)
head(wgcna_counts_filtered_long)## # A tibble: 6 × 3
## Gene Colony Counts
## <chr> <fct> <int>
## 1 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF13B 61
## 2 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF13D 73
## 3 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF14B 62
## 4 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF14C 51
## 5 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF15B 31
## 6 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF15D 37wgcna_counts_filtered_long <- wgcna_counts_filtered_long %>%
separate(Colony, into = c('Origin', 'Colony.number'), sep = 2)
head(wgcna_counts_filtered_long)## # A tibble: 6 × 4
## Gene Origin Colony.number Counts
## <chr> <chr> <chr> <int>
## 1 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF 13B 61
## 2 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF 13D 73
## 3 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF 14B 62
## 4 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF 14C 51
## 5 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF 15B 31
## 6 Pocillopora_acuta_HIv2___RNAseq.g27841.t1 RF 15D 37library(stringr)
wgcna_counts_filtered_long$Colony <- as.numeric(str_extract(wgcna_counts_filtered_long$Colony.number, "[0-9]+"))
wgcna_counts_filtered_long<-wgcna_counts_filtered_long %>%
mutate(Treatment = trimws(str_remove(wgcna_counts_filtered_long$Colony.number, "(\\s+[A-Za-z]+)?[0-9-]+")))
head(wgcna_counts_filtered_long)## # A tibble: 6 × 6
## Gene Origin Colony.number Counts Colony Treatment
## <chr> <chr> <chr> <int> <dbl> <chr>
## 1 Pocillopora_acuta_HIv2___RNAseq.… RF 13B 61 13 B
## 2 Pocillopora_acuta_HIv2___RNAseq.… RF 13D 73 13 D
## 3 Pocillopora_acuta_HIv2___RNAseq.… RF 14B 62 14 B
## 4 Pocillopora_acuta_HIv2___RNAseq.… RF 14C 51 14 C
## 5 Pocillopora_acuta_HIv2___RNAseq.… RF 15B 31 15 B
## 6 Pocillopora_acuta_HIv2___RNAseq.… RF 15D 37 15 Dwgcna_counts_filtered_long$Origin <- as.factor(wgcna_counts_filtered_long$Origin)
wgcna_counts_filtered_long$Treatment <- as.factor(wgcna_counts_filtered_long$Treatment)wgcna_counts_filtered_long <- wgcna_counts_filtered_long %>%
mutate(Treatment2 = ifelse(Treatment == "A" | Treatment == "B", "Variable",
ifelse(Treatment == "C" | Treatment == "D", "Stable", NA)))
wgcna_counts_filtered_long$Treatment2 <- as.factor(wgcna_counts_filtered_long$Treatment2)##SLC4A7
wgcna_counts_filtered_long_SLC4A7<- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g7402.t1")library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:IRanges':
##
## collapse## The following object is masked from 'package:dplyr':
##
## collapselibrary(emmeans)
SLC4A7.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_SLC4A7, na.action=na.exclude)
car::Anova(SLC4A7.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 11.9704 1 0.0005405 ***
## Origin 0.4869 1 0.4853006
## Treatment 1.0445 3 0.7904802
## Origin:Treatment 1.3516 3 0.7169263
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(SLC4A7.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 654 61.5 19 525 783
## RS 453 54.2 19 339 566
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 201 80.8 19 2.491 0.0221
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
SLC4A7_sum<-summarySE(wgcna_counts_filtered_long_SLC4A7, measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
SLC4A7_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 710.3636 265.8632 80.16078 178.6093
## 2 RF Variable 11 609.8182 241.3930 72.78272 162.1700
## 3 RS Stable 12 463.8333 248.6454 71.77773 157.9817
## 4 RS Variable 12 469.8333 166.4407 48.04730 105.7514Figure
pd<- position_dodge(0.2)
SLC4A7_fig<-ggplot(data=SLC4A7_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_SLC4A7,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(SLC4A7~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statSLC4A7_fig
##SLC4A3
wgcna_counts_filtered_long_SLC4A3 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___TS.g27873.t1")
wgcna_counts_filtered_long_SLC4A3 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 52 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 65 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 105 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 52 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 22 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 50 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 44 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 94 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 58 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 60 18 D Stable
## # ℹ 36 more rowsSLC4A3.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_SLC4A3 , na.action=na.exclude)
car::Anova(SLC4A3.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 24.6847 1 6.752e-07 ***
## Origin 1.7488 1 0.1860
## Treatment 0.4857 3 0.9220
## Origin:Treatment 1.1274 3 0.7705
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(SLC4A3.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 73.1 5.98 19 60.6 85.6
## RS 58.5 5.28 19 47.4 69.6
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 14.6 7.81 19 1.866 0.0776
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
SLC4A3_sum<-summarySE(wgcna_counts_filtered_long_SLC4A3 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
SLC4A3_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 70.63636 21.65767 6.530032 14.54982
## 2 RF Variable 11 71.72727 27.65896 8.339491 18.58154
## 3 RS Stable 12 61.50000 22.71763 6.558016 14.43410
## 4 RS Variable 12 56.33333 17.55166 5.066726 11.15179Figure
pd<- position_dodge(0.2)
SLC4A3_fig<-ggplot(data=SLC4A3_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_SLC4A3 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(SLC4A3 ~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statSLC4A3_fig
##NHE3
wgcna_counts_filtered_long_NHE3<- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g24868.t1")
wgcna_counts_filtered_long_NHE3## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 88 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 93 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 143 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 96 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 123 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 136 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 87 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 111 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 95 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 87 18 D Stable
## # ℹ 36 more rowsNHE3.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_NHE3, na.action=na.exclude)
car::Anova(NHE3.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 45.5717 1 1.472e-11 ***
## Origin 1.1334 1 0.2870
## Treatment 5.3488 3 0.1480
## Origin:Treatment 5.9678 3 0.1132
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(NHE3.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 111 6.90 19 96.9 126
## RS 102 6.19 19 89.0 115
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 9.42 8.58 19 1.098 0.2861
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
NHE3_sum<-summarySE(wgcna_counts_filtered_long_NHE3, measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
NHE3_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 105.18182 15.16455 4.572284 10.18768
## 2 RF Variable 11 121.63636 33.62818 10.139278 22.59172
## 3 RS Stable 12 106.66667 22.81281 6.585491 14.49457
## 4 RS Variable 12 97.83333 22.00757 6.353040 13.98295Figure
pd<- position_dodge(0.2)
NHE3_fig<-ggplot(data=NHE3_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_NHE3,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(NHE3~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statNHE3_fig
##CA1
wgcna_counts_filtered_long_CA1<- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___TS.g12304.t1")
wgcna_counts_filtered_long_CA1## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 2854 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 4635 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 2949 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 4681 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 7704 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 8665 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 3948 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 3887 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 6896 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 6597 18 D Stable
## # ℹ 36 more rowsCA1.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_CA1, na.action=na.exclude)
car::Anova(CA1.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 7.3722 1 0.006624 **
## Origin 0.7737 1 0.379079
## Treatment 2.6212 3 0.453784
## Origin:Treatment 2.3227 3 0.508194
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(CA1.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 5178 542 19 4045 6312
## RS 2568 479 19 1565 3571
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 2611 704 19 3.709 0.0015
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
CA1_sum<-summarySE(wgcna_counts_filtered_long_CA1, measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
CA1_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 5970.636 2621.084 790.2864 1760.8679
## 2 RF Variable 11 4733.455 1986.062 598.8204 1334.2549
## 3 RS Stable 12 2653.583 1847.177 533.2340 1173.6400
## 4 RS Variable 12 2775.250 1463.636 422.5154 929.9501Figure
pd<- position_dodge(0.2)
CA1_fig<-ggplot(data=CA1_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_CA1,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(CA1~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statCA1_fig
##CA2
wgcna_counts_filtered_long_CA2<- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g13824.t1")
wgcna_counts_filtered_long_CA2## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 139 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 164 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 169 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 249 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 174 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 232 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 232 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 204 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 244 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 338 18 D Stable
## # ℹ 36 more rowsCA2.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_CA2, na.action=na.exclude)
car::Anova(CA2.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 12.1927 1 0.0004798 ***
## Origin 8.6459 1 0.0032780 **
## Treatment 2.6667 3 0.4459141
## Origin:Treatment 1.7689 3 0.6217181
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(CA2.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 148.05 16.4 19 113.8 182.3
## RS -4.05 15.0 19 -35.5 27.4
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 152 18.2 19 8.346 <.0001
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
CA2_sum<-summarySE(wgcna_counts_filtered_long_CA2, measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
CA2_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 150.181818 98.679094 29.752866 66.293518
## 2 RF Variable 11 131.272727 70.816793 21.352066 47.575369
## 3 RS Stable 12 10.583333 9.894519 2.856302 6.286678
## 4 RS Variable 12 9.916667 8.106769 2.340223 5.150795Figure
pd<- position_dodge(0.2)
CA2_fig<-ggplot(data=CA2_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_CA2,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(CA2~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statCA2_fig
##VHA
wgcna_counts_filtered_long_VHA<- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g23064.t1")
wgcna_counts_filtered_long_VHA## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 28 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 19 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 26 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 21 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 35 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 30 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 30 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 21 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 16 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 17 18 D Stable
## # ℹ 36 more rowsVHA.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_VHA, na.action=na.exclude)
car::Anova(VHA.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 17.3800 1 3.06e-05 ***
## Origin 0.5371 1 0.4636
## Treatment 2.5674 3 0.4632
## Origin:Treatment 2.9805 3 0.3946
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(VHA.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 24.6 2.29 19 19.8 29.3
## RS 24.9 2.01 19 20.7 29.1
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -0.353 3.02 19 -0.117 0.9080
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
VHA_sum<-summarySE(wgcna_counts_filtered_long_VHA, measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
VHA_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 21.54545 5.317210 1.603199 3.572151
## 2 RF Variable 11 27.36364 5.427204 1.636364 3.646045
## 3 RS Stable 12 26.75000 13.784873 3.979350 8.758491
## 4 RS Variable 12 23.41667 7.025387 2.028054 4.463718Figure
pd<- position_dodge(0.2)
VHA_fig<-ggplot(data=VHA_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_VHA,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(VHA~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statVHA_fig
##HSP90
wgcna_counts_filtered_long_HSP90<- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g6656.t1")
wgcna_counts_filtered_long_HSP90## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 908 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 778 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 1301 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 891 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 1043 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 948 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 1547 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 906 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 1425 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 1568 18 D Stable
## # ℹ 36 more rowsHSP90.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_HSP90, na.action=na.exclude)
car::Anova(HSP90.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 34.6376 1 3.972e-09 ***
## Origin 0.5404 1 0.4623
## Treatment 2.0956 3 0.5528
## Origin:Treatment 0.7946 3 0.8508
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(HSP90.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 1345 92.2 19 1152 1538
## RS 1102 80.8 19 933 1271
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 243 123 19 1.985 0.0618
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
HSP90_sum<-summarySE(wgcna_counts_filtered_long_HSP90, measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
HSP90_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 1218.1818 382.7918 115.41607 257.1630
## 2 RF Variable 11 1426.5455 363.5721 109.62111 244.2511
## 3 RS Stable 12 969.3333 311.6707 89.97157 198.0261
## 4 RS Variable 12 1274.5833 397.4287 114.72777 252.5141Figure
pd<- position_dodge(0.2)
HSP90_fig<-ggplot(data=HSP90_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_HSP90,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(HSP90~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statHSP90_fig
##HIF1A
wgcna_counts_filtered_long_HIF1A <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g3039.t1")
wgcna_counts_filtered_long_HIF1A ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 172 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 158 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 99 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 159 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 127 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 164 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 93 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 128 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 104 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 213 18 D Stable
## # ℹ 36 more rowsHIF1A.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_HIF1A , na.action=na.exclude)
car::Anova(HIF1A.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 23.0105 1 1.611e-06 ***
## Origin 0.0263 1 0.8712
## Treatment 2.1843 3 0.5350
## Origin:Treatment 2.0693 3 0.5582
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(HIF1A.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 146 10.55 19 124 169
## RS 174 9.41 19 155 194
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -28 13.4 19 -2.097 0.0496
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
HIF1A_sum<-summarySE(wgcna_counts_filtered_long_HIF1A , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
HIF1A_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 151.5455 28.24310 8.515615 18.97397
## 2 RF Variable 11 141.2727 33.79672 10.190094 22.70494
## 3 RS Stable 12 192.0833 51.70627 14.926313 32.85259
## 4 RS Variable 12 165.1667 34.39565 9.929168 21.85395Figure
pd<- position_dodge(0.2)
HIF1A_fig<-ggplot(data=HIF1A_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_HIF1A ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(HIF1A ~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statHIF1A_fig
HSP70
wgcna_counts_filtered_long_HSP70 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g10659.t1")
wgcna_counts_filtered_long_HSP70 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 57 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 57 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 71 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 36 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 74 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 78 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 92 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 59 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 54 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 36 18 D Stable
## # ℹ 36 more rowsHSP70.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_HSP70 , na.action=na.exclude)
car::Anova(HSP70.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 10.4968 1 0.001196 **
## Origin 0.0701 1 0.791247
## Treatment 1.9098 3 0.591333
## Origin:Treatment 8.0987 3 0.044015 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(HSP70.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 63.1 10.78 19 40.5 85.6
## RS 103.0 9.72 19 82.7 123.3
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -40 13.1 19 -3.043 0.0067
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
HSP70_sum<-summarySE(wgcna_counts_filtered_long_HSP70 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
HSP70_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 59.09091 15.70003 4.733737 10.54742
## 2 RF Variable 11 71.63636 20.51474 6.185427 13.78199
## 3 RS Stable 12 116.00000 64.30750 18.563976 40.85904
## 4 RS Variable 12 95.33333 31.31463 9.039755 19.89637Figure
pd<- position_dodge(0.2)
HSP70_fig<-ggplot(data=HSP70_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_HSP70 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(HSP70 ~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statHSP70_fig
PRKCD
wgcna_counts_filtered_long_PRKCD <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g25259.t1")
wgcna_counts_filtered_long_PRKCD ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 7 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 9 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 47 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 18 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 25 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 62 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 22 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 35 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 28 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 23 18 D Stable
## # ℹ 36 more rowsPRKCD.lme <- lme(Counts~Origin*Treatment, random = ~1|Colony, data=wgcna_counts_filtered_long_PRKCD , na.action=na.exclude)
car::Anova(PRKCD.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 9.0491 1 0.002628 **
## Origin 3.7499 1 0.052810 .
## Treatment 5.4359 3 0.142525
## Origin:Treatment 2.4700 3 0.480731
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(PRKCD.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 25.40 3.25 19 18.60 32.2
## RS 7.18 2.91 19 1.09 13.3
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 18.2 4.04 19 4.508 0.0002
##
## Results are averaged over the levels of: Treatment
## Degrees-of-freedom method: containmentlibrary(Rmisc)
PRKCD_sum<-summarySE(wgcna_counts_filtered_long_PRKCD , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
PRKCD_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 25.000000 17.487138 5.272571 11.748019
## 2 RF Variable 11 26.636364 14.955084 4.509128 10.046962
## 3 RS Stable 12 7.333333 7.749878 2.237197 4.924037
## 4 RS Variable 12 7.083333 7.403419 2.137183 4.703908Figure
pd<- position_dodge(0.2)
PRKCD_fig<-ggplot(data=PRKCD_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(PRKCD ~expression))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statPRKCD_fig
compare_figs<-cowplot::plot_grid(SLC4A7_fig, NHE3_fig, CA1_fig, CA2_fig, HSP70_fig,HIF1A_fig, nrow=3)
compare_figs
Biomineralization toolkit present in modules
biomin <-read.csv("../../output/Biomin_blast_Pocillopora_acuta_best_hit.csv")wgcnamod <-read.csv("../../output/WGCNA/WGCNA_ModuleMembership.csv")
wgcnamod<- plyr::rename(wgcnamod, c("X"="Pocillopora_acuta_best_hit"))biomin_mod <- merge(biomin, wgcnamod, by=c("Pocillopora_acuta_best_hit"), all=F)
head(biomin_mod)## Pocillopora_acuta_best_hit accessionnumber.geneID
## 1 Pocillopora_acuta_HIv2___RNAseq.g10093.t2 XP_022804785.1
## 2 Pocillopora_acuta_HIv2___RNAseq.g11609.t1 P33_g8985
## 3 Pocillopora_acuta_HIv2___RNAseq.g13172.t1 JR972076.1
## 4 Pocillopora_acuta_HIv2___RNAseq.g13172.t1 Gene:g13552
## 5 Pocillopora_acuta_HIv2___RNAseq.g13172.t1 aug_v2a.06327.t1
## 6 Pocillopora_acuta_HIv2___RNAseq.g13823.t1 PFX18785.1
## definition
## 1 thioredoxin reductase 1, cytoplasmic-like [Stylophora pistillata]
## 2 Flagellar associated protein
## 3 Acidic skeletal organic matrix protein (Acidic SOMP)
## 4 Acidic SOMP (Full-Length p27)
## 5 SAARP3
## 6 Mucin-4 [Stylophora pistillata]
## Ref substanceBXH
## 1 Peled et al., 2020 Pocillopora_acuta_HIv2___RNAseq.g10093.t2
## 2 Drake et al., 2013 Pocillopora_acuta_HIv2___RNAseq.g11609.t1
## 3 Ramos-Silva et al., 2013 Pocillopora_acuta_HIv2___RNAseq.g13172.t1
## 4 Mummadisetti et al., 2021 Pocillopora_acuta_HIv2___RNAseq.g13172.t1
## 5 Takeuchi et al., 2016 Pocillopora_acuta_HIv2___RNAseq.g13172.t1
## 6 Peled et al., 2020 Pocillopora_acuta_HIv2___RNAseq.g13823.t1
## geneSymbol moduleColor
## 1 Pocillopora_acuta_HIv2___Sc0000021 brown
## 2 Pocillopora_acuta_HIv2___Sc0000013 turquoise
## 3 Pocillopora_acuta_HIv2___Sc0000004 red
## 4 Pocillopora_acuta_HIv2___Sc0000004 red
## 5 Pocillopora_acuta_HIv2___Sc0000004 red
## 6 Pocillopora_acuta_HIv2___Sc0000005 pink
## GO.terms
## 1 GO:0000003,GO:0000302,GO:0000305,GO:0001650,GO:0001704,GO:0001707,GO:0001887,GO:0001890,GO:0003006,GO:0003674,GO:0003824,GO:0004791,GO:0005488,GO:0005515,GO:0005575,GO:0005622,GO:0005623,GO:0005634,GO:0005654,GO:0005730,GO:0005737,GO:0005739,GO:0005783,GO:0005829,GO:0006082,GO:0006139,GO:0006518,GO:0006520,GO:0006575,GO:0006725,GO:0006732,GO:0006733,GO:0006739,GO:0006749,GO:0006753,GO:0006790,GO:0006793,GO:0006796,GO:0006807,GO:0006950,GO:0006979,GO:0007154,GO:0007165,GO:0007275,GO:0007369,GO:0007498,GO:0008150,GO:0008152,GO:0008283,GO:0009056,GO:0009069,GO:0009117,GO:0009611,GO:0009628,GO:0009636,GO:0009653,GO:0009790,GO:0009888,GO:0009987,GO:0010035,GO:0010038,GO:0010269,GO:0010941,GO:0010942,GO:0012505,GO:0015036,GO:0015949,GO:0016043,GO:0016174,GO:0016209,GO:0016259,GO:0016491,GO:0016651,GO:0016667,GO:0016668,GO:0016999,GO:0017001,GO:0017144,GO:0018996,GO:0019216,GO:0019222,GO:0019362,GO:0019637,GO:0019725,GO:0019752,GO:0022404,GO:0022414,GO:0022607,GO:0023052,GO:0031974,GO:0031981,GO:0032501,GO:0032502,GO:0033554,GO:0033797,GO:0034599,GO:0034641,GO:0036295,GO:0036296,GO:0036477,GO:0042221,GO:0042303,GO:0042395,GO:0042493,GO:0042537,GO:0042592,GO:0042737,GO:0042743,GO:0042744,GO:0042802,GO:0042803,GO:0043025,GO:0043167,GO:0043169,GO:0043226,GO:0043227,GO:0043228,GO:0043229,GO:0043231,GO:0043232,GO:0043233,GO:0043436,GO:0043603,GO:0043933,GO:0044085,GO:0044237,GO:0044238,GO:0044248,GO:0044281,GO:0044297,GO:0044422,GO:0044424,GO:0044428,GO:0044444,GO:0044446,GO:0044452,GO:0044464,GO:0045340,GO:0045454,GO:0046483,GO:0046496,GO:0046688,GO:0046872,GO:0046914,GO:0046983,GO:0048332,GO:0048513,GO:0048518,GO:0048522,GO:0048598,GO:0048608,GO:0048646,GO:0048678,GO:0048729,GO:0048731,GO:0048856,GO:0050664,GO:0050789,GO:0050794,GO:0050896,GO:0051186,GO:0051187,GO:0051259,GO:0051262,GO:0051716,GO:0055086,GO:0055093,GO:0055114,GO:0061458,GO:0065003,GO:0065007,GO:0065008,GO:0070013,GO:0070276,GO:0070482,GO:0070887,GO:0070995,GO:0071241,GO:0071248,GO:0071280,GO:0071453,GO:0071455,GO:0071704,GO:0071840,GO:0072524,GO:0072593,GO:0080090,GO:0097237,GO:0097458,GO:0098623,GO:0098625,GO:0098626,GO:0098754,GO:0098869,GO:1901360,GO:1901564,GO:1901605,GO:1901700,GO:1990748
## 2 -
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## GO.description GS.Flat GS.Slope p.GS.Flat
## 1 thioredoxin-disulfide reductase activity 0.57178848 -0.57178848 3.311055e-05
## 2 - -0.29586493 0.29586493 4.589336e-02
## 3 <NA> 0.35628512 -0.35628512 1.508700e-02
## 4 <NA> 0.35628512 -0.35628512 1.508700e-02
## 5 <NA> 0.35628512 -0.35628512 1.508700e-02
## 6 <NA> -0.05455251 0.05455251 7.187880e-01
## p.GS.Slope A.brown p.A.brown A.magenta p.A.magenta A.red
## 1 3.311055e-05 0.7005073 5.973619e-08 -0.3738439 1.048844e-02 0.2901298
## 2 4.589336e-02 -0.4291375 2.921081e-03 0.3115539 3.505853e-02 -0.3452015
## 3 1.508700e-02 0.4914202 5.241968e-04 -0.6288605 2.864308e-06 0.6673892
## 4 1.508700e-02 0.4914202 5.241968e-04 -0.6288605 2.864308e-06 0.6673892
## 5 1.508700e-02 0.4914202 5.241968e-04 -0.6288605 2.864308e-06 0.6673892
## 6 7.187880e-01 0.0972208 5.203783e-01 -0.3252127 2.743096e-02 0.3709019
## p.A.red A.turquoise p.A.turquoise A.purple p.A.purple A.green
## 1 5.047695e-02 -0.43323233 2.634293e-03 0.6984202 6.792759e-08 0.4574538
## 2 1.879503e-02 0.58815287 1.720729e-05 -0.1784560 2.353887e-01 -0.1306835
## 3 4.071016e-07 -0.13892006 3.571825e-01 0.1198762 4.274677e-01 0.2378899
## 4 4.071016e-07 -0.13892006 3.571825e-01 0.1198762 4.274677e-01 0.2378899
## 5 4.071016e-07 -0.13892006 3.571825e-01 0.1198762 4.274677e-01 0.2378899
## 6 1.116224e-02 0.08164806 5.895928e-01 -0.1391597 3.563456e-01 0.1614616
## p.A.green A.lightcyan p.A.lightcyan A.pink p.A.pink A.blue
## 1 0.001391986 -0.3508191 1.682948e-02 0.1707384 2.565893e-01 0.12358439
## 2 0.386672688 0.1196505 4.283449e-01 -0.1522331 3.125037e-01 -0.58598406
## 3 0.111386103 -0.6473842 1.159989e-06 0.7188738 1.835918e-08 0.07448551
## 4 0.111386103 -0.6473842 1.159989e-06 0.7188738 1.835918e-08 0.07448551
## 5 0.111386103 -0.6473842 1.159989e-06 0.7188738 1.835918e-08 0.07448551
## 6 0.283719167 -0.5276145 1.646022e-04 0.6417477 1.537006e-06 -0.02286640
## p.A.blue A.salmon p.A.salmon A.midnightblue p.A.midnightblue
## 1 4.132051e-01 0.1178467 0.435389343 0.2439890 0.10224333
## 2 1.880492e-05 0.1907995 0.204028320 0.2258109 0.13131383
## 3 6.227429e-01 0.4254256 0.003204458 0.2691022 0.07053914
## 4 6.227429e-01 0.4254256 0.003204458 0.2691022 0.07053914
## 5 6.227429e-01 0.4254256 0.003204458 0.2691022 0.07053914
## 6 8.801027e-01 0.2940377 0.047315397 0.2906592 0.05003895
## A.black p.A.black A.cyan p.A.cyan A.yellow p.A.yellow
## 1 -0.28430645 0.05550307 0.04904562 0.7461773 0.05522073 0.7154873547
## 2 -0.18825739 0.21023361 0.07386502 0.6256510 -0.14392338 0.3399558326
## 3 0.09618758 0.52484209 0.16699226 0.2673276 -0.38010677 0.0091694889
## 4 0.09618758 0.52484209 0.16699226 0.2673276 -0.38010677 0.0091694889
## 5 0.09618758 0.52484209 0.16699226 0.2673276 -0.38010677 0.0091694889
## 6 0.02103556 0.88964060 -0.13338389 0.3768501 -0.46998526 0.0009821983
## A.tan p.A.tan
## 1 0.2648346 0.075293267
## 2 0.2613466 0.079363532
## 3 -0.2055446 0.170565420
## 4 -0.2055446 0.170565420
## 5 -0.2055446 0.170565420
## 6 -0.3805358 0.009084622plyr::count(biomin_mod, "moduleColor")## moduleColor freq
## 1 black 3
## 2 blue 36
## 3 brown 17
## 4 cyan 2
## 5 green 3
## 6 magenta 1
## 7 pink 7
## 8 red 18
## 9 salmon 6
## 10 tan 3
## 11 turquoise 21
## 12 yellow 10Format GO terms to remove dashes and quotes and separate by semicolons (replace , with ;) in GOs column
biomin_mod$GO.terms <- gsub(",", ";", biomin_mod$GO.terms)
biomin_mod$GO.terms <- gsub('"', "", biomin_mod$GO.terms)
biomin_mod$GO.terms <- gsub("-", NA, biomin_mod$GO.terms)GO terms
### Generate vector with names of all genes
ALL.vector <- c(geneInfo$gene_id)### Generate length vector for all genes
LENGTH.vector <- as.integer(geneInfo$Length)ID.vector_biomin <- biomin_mod %>%
#filter(moduleColor=="black")%>%
#get_rows(.data[[module]]))%>%
pull(Pocillopora_acuta_best_hit)
##Get a list of GO Terms for each module
GO.terms_biomin <- biomin_mod %>%
#filter(moduleColor=="black")%>%
#filter(get_rows(.data[[module]]))%>%
dplyr::select(GO.terms,Pocillopora_acuta_best_hit)##Format to have one goterm per row with gene ID repeated
split <- strsplit(as.character(GO.terms_biomin$GO.terms), ";")
split2 <- data.frame(v1 = rep.int(GO.terms_biomin$Pocillopora_acuta_best_hit, sapply(split, length)), v2 = unlist(split))
colnames(split2) <- c("Pocillopora_acuta_best_hit", "GO.terms")
GO.terms_biomin<-split2
head(GO.terms_biomin)## Pocillopora_acuta_best_hit GO.terms
## 1 Pocillopora_acuta_HIv2___RNAseq.g10093.t2 GO:0000003
## 2 Pocillopora_acuta_HIv2___RNAseq.g10093.t2 GO:0000302
## 3 Pocillopora_acuta_HIv2___RNAseq.g10093.t2 GO:0000305
## 4 Pocillopora_acuta_HIv2___RNAseq.g10093.t2 GO:0001650
## 5 Pocillopora_acuta_HIv2___RNAseq.g10093.t2 GO:0001704
## 6 Pocillopora_acuta_HIv2___RNAseq.g10093.t2 GO:0001707#GO.terms_biomin_sub <- GO.terms_biomin%>%
#filter(Pocillopora_acuta_best_hit==c("Pocillopora_acuta_HIv2___RNAseq.g15280.t1","Pocillopora_acuta_HIv2___RNAseq.g7402.t1"))
#GO.terms_biomin_sub##Construct list of genes with 1 for genes in module and 0 for genes not in the module
gene.vector=as.integer(ALL.vector %in% ID.vector)
names(gene.vector)<-ALL.vector#set names
#weight gene vector by bias for length of gene
pwf<-nullp(gene.vector, ID.vector, bias.data=LENGTH.vector) ## Warning in pcls(G): initial point very close to some inequality constraints
#run goseq using Wallenius method for all categories of GO terms
GO.wall<-goseq(pwf, ID.vector, gene2cat=GO.terms_biomin, test.cats=c("GO:BP", "GO:MF", "GO:CC"), method="Wallenius", use_genes_without_cat=TRUE)## Using manually entered categories.## Calculating the p-values...## 'select()' returned 1:1 mapping between keys and columnsGO <- GO.wall[order(GO.wall$over_represented_pvalue),]
colnames(GO)[1] <- "GOterm"#adjust p-values
GO$bh_adjust <- p.adjust(GO$over_represented_pvalue, method="BH") #add adjusted p-values#Filtering for p < 0.01
GO <- GO %>%
#dplyr::filter(bh_adjust<0.05) %>%
dplyr::arrange(., ontology, bh_adjust)
#Write file of results
write.csv(GO, file = "../../output/WGCNA/GO_analysis/goseq_pattern_biomin.csv")#add vector for terms of interest to reduce number of GO terms - NOT using this to look at individual modules for exploratory purposes
keywords<-c("metabolism", "carbon","bicarbonate", "apoptosis", "death", "symbiosis", "regulation of cell communication", "trans membrane transport", "transmembrane", "organic substance transport", "inorganic substance transport","response to stress", "antioxidant", "calcification","biomineralization", "heat","transporters","proton transport","ion transport","acid-base", "homeostasis")go_results <-read.csv("../../output/WGCNA/GO_analysis/goseq_pattern_biomin.csv") go_results<-go_results%>%
filter(ontology=="BP")%>%
filter(bh_adjust != "NA") %>%
#filter(numInCat>10)%>%
arrange(., bh_adjust)
head(go_results)## X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 GO:0006325 0.009864096 1 2
## 2 2 GO:0016570 0.009864096 1 2
## 3 3 GO:0051276 0.009864096 1 2
## 4 4 GO:0060537 0.012664490 1 2
## 5 5 GO:0000278 0.012973357 1 2
## 6 6 GO:0007049 0.012973357 1 2
## numInCat term ontology bh_adjust
## 1 2 chromatin organization BP 0.8862339
## 2 2 histone modification BP 0.8862339
## 3 2 chromosome organization BP 0.8862339
## 4 2 muscle tissue development BP 0.8862339
## 5 2 mitotic cell cycle BP 0.8862339
## 6 2 cell cycle BP 0.8862339library(rrvgo)
#Reduce/collapse GO term set with the rrvgo package
simMatrix <- calculateSimMatrix(go_results$GOterm,
orgdb="org.Ce.eg.db", #c. elegans database
ont="BP",
method="Rel")## preparing gene to GO mapping data...## preparing IC data...#calculate similarity
scores <- setNames(-log(go_results$bh_adjust), go_results$GOterm)
reducedTerms <- reduceSimMatrix(simMatrix,
scores,
threshold=0.7,
orgdb="org.Ce.eg.db")
head(reducedTerms)## go cluster parent score size
## GO:0030097 GO:0030097 30 GO:0030097 0.1207744 0
## GO:0002376 GO:0002376 25 GO:0002376 0.1207744 24
## GO:0032501 GO:0032501 19 GO:0032501 0.1207744 2
## GO:0044260 GO:0044260 7 GO:0044260 0.1207744 12
## GO:0006338 GO:0006338 1 GO:0006338 0.1207744 59
## GO:0008210 GO:0008210 26 GO:0008210 0.1207744 2
## term
## GO:0030097 hemopoiesis
## GO:0002376 immune system process
## GO:0032501 multicellular organismal process
## GO:0044260 cellular macromolecule metabolic process
## GO:0006338 chromatin remodeling
## GO:0008210 estrogen metabolic process
## parentTerm termUniqueness
## GO:0030097 hemopoiesis 0.9252500
## GO:0002376 immune system process 1.0000000
## GO:0032501 multicellular organismal process 0.9398698
## GO:0044260 cellular macromolecule metabolic process 0.9288776
## GO:0006338 chromatin remodeling 0.9418763
## GO:0008210 estrogen metabolic process 0.9288750
## termUniquenessWithinCluster termDispensability
## GO:0030097 1.0000000 0
## GO:0002376 1.0000000 0
## GO:0032501 0.6136250 0
## GO:0044260 0.6100000 0
## GO:0006338 0.5873333 0
## GO:0008210 0.5872500 0#keep only the goterms from the reduced list
go_results<-go_results%>%
filter(GOterm %in% reducedTerms$go)
#add in parent terms to list of go terms
go_results$ParentTerm<-reducedTerms$parentTerm[match(go_results$GOterm, reducedTerms$go)]
go_results<-go_results %>%
mutate(Factor = "Biomin")
head(go_results)## X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 GO:0006325 0.009864096 1 2
## 2 2 GO:0016570 0.009864096 1 2
## 3 3 GO:0051276 0.009864096 1 2
## 4 4 GO:0060537 0.012664490 1 2
## 5 5 GO:0000278 0.012973357 1 2
## 6 6 GO:0007049 0.012973357 1 2
## numInCat term ontology bh_adjust
## 1 2 chromatin organization BP 0.8862339
## 2 2 histone modification BP 0.8862339
## 3 2 chromosome organization BP 0.8862339
## 4 2 muscle tissue development BP 0.8862339
## 5 2 mitotic cell cycle BP 0.8862339
## 6 2 cell cycle BP 0.8862339
## ParentTerm Factor
## 1 chromatin remodeling Biomin
## 2 macromolecule deacylation Biomin
## 3 chromatin remodeling Biomin
## 4 muscle tissue development Biomin
## 5 microtubule cytoskeleton organization involved in mitosis Biomin
## 6 microtubule cytoskeleton organization involved in mitosis Biominwrite.csv(go_results, "../../output/Biomineralization_goterms.csv")go_results<-go_results%>%
filter(grepl(paste(keywords, collapse="|"), ParentTerm))#plot significantly enriched GO terms by Slim Category faceted by slim term
GO.plot_biomin <- ggplot(go_results, aes(x = ontology, y = term)) +
geom_point(aes(size=bh_adjust)) +
scale_size(name="Over rep. p-value", trans="reverse", range=c(1,3))+
facet_grid(ParentTerm ~ ., scales = "free", labeller = label_wrap_gen(width = 5, multi_line = TRUE))+
theme_bw() + theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"),
strip.text.y = element_text(angle=0, size = 10),
strip.text.x = element_text(size = 20),
axis.text = element_text(size = 8),
axis.title.x = element_blank(),
axis.title.y = element_blank())
GO.plot_biomin
***thioredoxin reductase 1
wgcna_counts_filtered_long_g10093 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g10093.t2")
wgcna_counts_filtered_long_g10093 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 59 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 69 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 86 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 75 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 49 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 63 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 80 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 57 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 98 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 73 18 D Stable
## # ℹ 36 more rowsg10093.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g10093 , na.action=na.exclude)
car::Anova(g10093.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 215.5517 1 < 2.2e-16 ***
## Origin 7.8998 1 0.004944 **
## Treatment2 1.0425 1 0.307243
## Origin:Treatment2 2.7791 1 0.095503 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g10093.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 66.5 3.57 19 59.0 73.9
## RS 44.8 3.43 19 37.6 52.0
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 21.6 4.47 23 4.836 0.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g10093_sum<-summarySE(wgcna_counts_filtered_long_g10093 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g10093_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 63.09091 12.70004 3.829205 8.532000
## 2 RF Variable 11 68.18182 19.91893 6.005782 13.381717
## 3 RS Stable 12 48.50000 12.57342 3.629634 7.988770
## 4 RS Variable 12 42.08333 12.44960 3.593889 7.910096Figure
pd<- position_dodge(0.2)
g10093_fig<-ggplot(data=g10093_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g10093,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Thioredoxin~reductase~1~expression))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg10093_fig
***Flagellar associated protein
wgcna_counts_filtered_long_g11609 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g11609.t1")
wgcna_counts_filtered_long_g11609 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 270 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 210 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 227 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 162 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 252 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 234 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 252 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 143 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 164 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 209 18 D Stable
## # ℹ 36 more rowsg11609.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g11609 , na.action=na.exclude)
car::Anova(g11609.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 55.8226 1 7.932e-14 ***
## Origin 10.3843 1 0.001271 **
## Treatment2 2.7386 1 0.097951 .
## Origin:Treatment2 9.2919 1 0.002302 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g11609.lme, list(pairwise ~ Treatment2:Origin), adjust = "tukey")
tukey3## $`emmeans of Treatment2, Origin`
## Treatment2 Origin emmean SE df lower.CL upper.CL
## Stable RF 213 28.5 19 153 273
## Variable RF 274 28.5 19 214 334
## Stable RS 337 27.3 19 280 395
## Variable RS 243 27.3 19 186 300
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Treatment2, Origin`
## 1 estimate SE df t.ratio p.value
## Stable RF - Variable RF -60.9 36.8 23 -1.655 0.3694
## Stable RF - Stable RS -124.1 38.5 23 -3.222 0.0184
## Stable RF - Variable RS -29.7 38.5 23 -0.772 0.8664
## Variable RF - Stable RS -63.2 38.5 23 -1.641 0.3763
## Variable RF - Variable RS 31.2 38.5 23 0.809 0.8494
## Stable RS - Variable RS 94.4 35.2 23 2.679 0.0601
##
## Degrees-of-freedom method: containment
## P value adjustment: tukey method for comparing a family of 4 estimateslibrary(Rmisc)
g11609_sum<-summarySE(wgcna_counts_filtered_long_g11609 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g11609_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 216.3636 68.28803 20.58961 45.87652
## 2 RF Variable 11 277.2727 116.48699 35.12215 78.25702
## 3 RS Stable 12 335.0833 113.94772 32.89387 72.39893
## 4 RS Variable 12 240.6667 69.27985 19.99937 44.01831Figure
pd<- position_dodge(0.2)
g11609_fig<-ggplot(data=g11609_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g11609 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Flagellar~associated~protein~expression))+
ggtitle(~turquoise)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg11609_fig
***Actin
wgcna_counts_filtered_long_g14505 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g14505.t1")
wgcna_counts_filtered_long_g14505 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 1418 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 1504 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 1744 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 1503 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 1805 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 1848 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 1405 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 1124 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 1361 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 1045 18 D Stable
## # ℹ 36 more rowsg14505.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g14505 , na.action=na.exclude)
car::Anova(g14505.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 134.3187 1 < 2.2e-16 ***
## Origin 7.2528 1 0.007079 **
## Treatment2 0.4459 1 0.504290
## Origin:Treatment2 2.6857 1 0.101254
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g14505.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 1631 95.6 19 1431 1832
## RS 1919 91.6 19 1727 2110
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -287 132 23 -2.170 0.0406
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g14505_sum<-summarySE(wgcna_counts_filtered_long_g14505 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g14505_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 1567.545 325.2954 98.08024 218.5364
## 2 RF Variable 11 1695.273 508.7410 153.39119 341.7769
## 3 RS Stable 12 2071.833 491.8505 141.98500 312.5069
## 4 RS Variable 12 1765.583 441.5165 127.45483 280.5262Figure
pd<- position_dodge(0.2)
g14505_fig<-ggplot(data=g14505_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g14505 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Actin~expression))+
ggtitle(~turquoise)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg14505_fig
***CA2
wgcna_counts_filtered_long_CA2<- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g13824.t1")
wgcna_counts_filtered_long_CA2## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 139 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 164 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 169 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 249 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 174 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 232 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 232 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 204 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 244 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 338 18 D Stable
## # ℹ 36 more rowsCA2.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_CA2, na.action=na.exclude)
car::Anova(CA2.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 96.0895 1 <2e-16 ***
## Origin 100.2939 1 <2e-16 ***
## Treatment2 2.1871 1 0.1392
## Origin:Treatment2 1.0620 1 0.3027
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(CA2.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 155.76 15.6 19 123.1 188
## RS -8.88 15.2 19 -40.8 23
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 165 14.9 23 11.032 <.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
CA2_sum<-summarySE(wgcna_counts_filtered_long_CA2, measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
CA2_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 150.181818 98.679094 29.752866 66.293518
## 2 RF Variable 11 131.272727 70.816793 21.352066 47.575369
## 3 RS Stable 12 10.583333 9.894519 2.856302 6.286678
## 4 RS Variable 12 9.916667 8.106769 2.340223 5.150795Figure
pd<- position_dodge(0.2)
CA2_fig<-ggplot(data=CA2_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_CA2,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(CA2~expression))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statCA2_fig
Poly [ADP-ribose] polymerase 11
wgcna_counts_filtered_long_g14663 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g14663.t1a")
wgcna_counts_filtered_long_g14663 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 61 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 54 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 19 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 24 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 23 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 49 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 17 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 13 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 16 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 10 18 D Stable
## # ℹ 36 more rowsg14663.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g14663 , na.action=na.exclude)
car::Anova(g14663.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 14.3519 1 0.0001516 ***
## Origin 0.5329 1 0.4653839
## Treatment2 0.0388 1 0.8437699
## Origin:Treatment2 0.0181 1 0.8928700
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g14663.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 25.5 5.73 19 13.5 37.5
## RS 32.8 5.52 19 21.3 44.4
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -7.28 7.13 23 -1.020 0.3184
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g14663_sum<-summarySE(wgcna_counts_filtered_long_g14663 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g14663_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 28.09091 15.23453 4.593384 10.234697
## 2 RF Variable 11 26.54545 14.34129 4.324063 9.634613
## 3 RS Stable 12 33.75000 28.81958 8.319496 18.311087
## 4 RS Variable 12 33.66667 30.09782 8.688492 19.123243Figure
pd<- position_dodge(0.2)
g14663_fig<-ggplot(data=g14663_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Poly~"[ADP-ribose]"~polymerase~11~expression))+
ggtitle(~turquoise)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg14663_fig
***solute carrier family 4 member gamma
wgcna_counts_filtered_long_g15280 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g15280.t1")
wgcna_counts_filtered_long_g15280 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 30 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 34 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 29 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 30 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 94 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 76 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 49 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 39 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 64 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 64 18 D Stable
## # ℹ 36 more rowsg15280.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g15280 , na.action=na.exclude)
car::Anova(g15280.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 86.8426 1 < 2.2e-16 ***
## Origin 8.6455 1 0.003279 **
## Treatment2 0.8180 1 0.365763
## Origin:Treatment2 0.8151 1 0.366627
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g15280.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 51.8 4.22 19 43.0 60.7
## RS 32.8 4.04 19 24.3 41.3
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 19 5.84 23 3.255 0.0035
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g15280_sum<-summarySE(wgcna_counts_filtered_long_g15280 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g15280_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 55.63636 24.46333 7.375972 16.43469
## 2 RF Variable 11 48.00000 20.63008 6.220202 13.85947
## 3 RS Stable 12 31.33333 16.77299 4.841946 10.65705
## 4 RS Variable 12 34.25000 16.87454 4.871259 10.72157Figure
pd<- position_dodge(0.2)
g15280_fig<-ggplot(data=g15280_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g15280 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Solute~carrier~family~4~member~gamma~expression))+
ggtitle(~pink)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg15280_fig
MAGUK p55 subfamily member 7-like
wgcna_counts_filtered_long_g15517 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g15517.t1")
wgcna_counts_filtered_long_g15517 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 0 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 14 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 27 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 15 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 12 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 14 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 15 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 0 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 9 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 29 18 D Stable
## # ℹ 36 more rowsg15517.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g15517 , na.action=na.exclude)
car::Anova(g15517.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 27.9449 1 1.248e-07 ***
## Origin 0.2433 1 0.6219
## Treatment2 0.2642 1 0.6073
## Origin:Treatment2 1.3103 1 0.2523
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g15517.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 13.5 1.95 19 9.47 17.6
## RS 14.8 1.86 19 10.85 18.6
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -1.2 2.69 23 -0.447 0.6589
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g15517_sum<-summarySE(wgcna_counts_filtered_long_g15517 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g15517_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 14.54545 9.490665 2.861543 6.375916
## 2 RF Variable 11 12.54545 10.103105 3.046201 6.787358
## 3 RS Stable 12 12.66667 7.992421 2.307213 5.078142
## 4 RS Variable 12 16.83333 8.912028 2.572681 5.662432Figure
pd<- position_dodge(0.2)
g15517_fig<-ggplot(data=g15517_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(MAGUK~p55~subfamily~member~7-like~expression))+
ggtitle(~cyan)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg15517_fig
CARP1
wgcna_counts_filtered_long_g16280 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g16280.t1")
wgcna_counts_filtered_long_g16280 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 103 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 165 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 384 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 232 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 309 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 396 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 225 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 228 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 187 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 112 18 D Stable
## # ℹ 36 more rowsg16280.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g16280 , na.action=na.exclude)
car::Anova(g16280.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 68.2145 1 <2e-16 ***
## Origin 0.1913 1 0.6619
## Treatment2 0.0758 1 0.7831
## Origin:Treatment2 0.0179 1 0.8934
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g16280.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 276 26.3 19 221 331
## RS 260 25.2 19 207 313
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 16.1 34.7 23 0.463 0.6477
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g16280_sum<-summarySE(wgcna_counts_filtered_long_g16280 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g16280_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 287.2727 124.33752 37.48917 83.53108
## 2 RF Variable 11 275.2727 90.98691 27.43359 61.12584
## 3 RS Stable 12 259.7500 119.95766 34.62879 76.21746
## 4 RS Variable 12 255.8333 115.64824 33.38477 73.47939Figure
pd<- position_dodge(0.2)
g16280_fig<-ggplot(data=g16280_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(CARP1~expression))+
ggtitle(~tan)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg16280_fig
PHD finger protein 21A-like
wgcna_counts_filtered_long_g1634 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g1634.t1")
wgcna_counts_filtered_long_g1634 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 65 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 80 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 90 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 79 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 95 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 95 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 76 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 78 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 91 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 79 18 D Stable
## # ℹ 36 more rowsg1634.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g1634 , na.action=na.exclude)
car::Anova(g1634.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 177.4960 1 <2e-16 ***
## Origin 0.5303 1 0.4665
## Treatment2 0.4221 1 0.5159
## Origin:Treatment2 0.0042 1 0.9480
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g1634.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 90.0 4.62 19 80.3 99.7
## RS 96.2 4.42 19 86.9 105.4
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -6.17 6.39 23 -0.965 0.3448
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g1634_sum<-summarySE(wgcna_counts_filtered_long_g1634 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g1634_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 87.00000 11.13553 3.357488 7.48095
## 2 RF Variable 11 93.00000 17.89413 5.395284 12.02144
## 3 RS Stable 12 93.58333 25.61412 7.394161 16.27444
## 4 RS Variable 12 98.75000 27.03911 7.805520 17.17983Figure
pd<- position_dodge(0.2)
g1634_fig<-ggplot(data=g1634_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(PHD~finger~protein~"21A"~like~expression))+
ggtitle(~turquoise)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg1634_fig
digestive cysteine proteinase 1-like
wgcna_counts_filtered_long_g18103 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g18103.t1")
wgcna_counts_filtered_long_g18103 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 0 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 9 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 0 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 12 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 6 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 0 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 19 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 0 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 10 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 10 18 D Stable
## # ℹ 36 more rowsg18103.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g18103 , na.action=na.exclude)
car::Anova(g18103.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 14.6577 1 0.0001289 ***
## Origin 1.8257 1 0.1766413
## Treatment2 0.8843 1 0.3470176
## Origin:Treatment2 0.0057 1 0.9398383
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g18103.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 10.14 1.60 19 6.80 13.47
## RS 5.75 1.53 19 2.55 8.95
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 4.39 2.21 23 1.986 0.0590
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g18103_sum<-summarySE(wgcna_counts_filtered_long_g18103 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g18103_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 8.636364 6.313046 1.903455 4.241162
## 2 RF Variable 11 11.636364 10.032674 3.024965 6.740042
## 3 RS Stable 12 4.416667 7.216878 2.083333 4.585386
## 4 RS Variable 12 7.083333 5.822501 1.680811 3.699440Figure
pd<- position_dodge(0.2)
g18103_fig<-ggplot(data=g18103_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g18103 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(digestive~cysteine~proteinase~1-like~expression))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg18103_fig
endothelin-converting enzyme 1-like isoform X2
wgcna_counts_filtered_long_g19211 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g19211.t1")
wgcna_counts_filtered_long_g19211 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 46 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 56 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 51 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 30 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 43 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 27 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 37 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 47 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 54 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 67 18 D Stable
## # ℹ 36 more rowsg19211.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g19211 , na.action=na.exclude)
car::Anova(g19211.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 87.3886 1 < 2e-16 ***
## Origin 3.4868 1 0.06186 .
## Treatment2 0.1819 1 0.66975
## Origin:Treatment2 0.8990 1 0.34305
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g19211.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 48.7 3.93 19 40.5 56.9
## RS 57.2 3.77 19 49.3 65.1
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -8.49 5.15 23 -1.647 0.1131
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g19211_sum<-summarySE(wgcna_counts_filtered_long_g19211 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g19211_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 48.18182 12.97550 3.912261 8.717060
## 2 RF Variable 11 50.90909 12.33251 3.718393 8.285096
## 3 RS Stable 12 59.16667 22.09415 6.378032 14.037954
## 4 RS Variable 12 53.50000 17.65065 5.095304 11.214688Figure
pd<- position_dodge(0.2)
g19211_fig<-ggplot(data=g19211_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(endothelin-converting~enzyme~1-like~isoform~X2~expression))+
ggtitle(~green)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg19211_fig
chymotrypsin-like elastase family member 1
wgcna_counts_filtered_long_g19288 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g19288.t1")
wgcna_counts_filtered_long_g19288 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 0 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 0 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 0 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 0 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 0 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 48 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 0 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 25 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 29 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 52 18 D Stable
## # ℹ 36 more rowsg19288.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g19288 , na.action=na.exclude)
car::Anova(g19288.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 6.7903 1 0.009165 **
## Origin 0.0148 1 0.903331
## Treatment2 2.4178 1 0.119960
## Origin:Treatment2 0.0599 1 0.806669
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g19288.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 9.58 4.58 19 -0.0019 19.2
## RS 11.54 4.40 19 2.3230 20.8
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -1.96 5.73 23 -0.342 0.7357
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g19288_sum<-summarySE(wgcna_counts_filtered_long_g19288 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g19288_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 13.272727 20.337605 6.132019 13.662989
## 2 RF Variable 11 3.363636 8.834848 2.663807 5.935332
## 3 RS Stable 12 15.833333 25.785244 7.443559 16.383162
## 4 RS Variable 12 8.083333 15.264089 4.406363 9.698340Figure
pd<- position_dodge(0.2)
g19288_fig<-ggplot(data=g19288_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(chymotrypsin-like~elastase~family~member~1~expression))+
ggtitle(~cyan)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg19288_fig
CUB and peptidase domain-containing protein 2-like
wgcna_counts_filtered_long_g21338 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g21338.t1")
wgcna_counts_filtered_long_g21338 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 6 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 0 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 7 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 9 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 25 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 20 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 7 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 0 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 15 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 32 18 D Stable
## # ℹ 36 more rowsg21338.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g21338 , na.action=na.exclude)
car::Anova(g21338.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 19.7037 1 9.043e-06 ***
## Origin 0.8360 1 0.3605
## Treatment2 0.0175 1 0.8947
## Origin:Treatment2 2.2965 1 0.1297
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g21338.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 14.3 2.66 19 8.71 19.8
## RS 14.3 2.56 19 8.99 19.7
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -0.0721 3.24 23 -0.022 0.9824
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g21338_sum<-summarySE(wgcna_counts_filtered_long_g21338 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g21338_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 13.00000 11.063453 3.335757 7.432529
## 2 RF Variable 11 13.45455 7.353416 2.217138 4.940092
## 3 RS Stable 12 17.83333 13.354150 3.855011 8.484822
## 4 RS Variable 12 11.08333 10.422514 3.008720 6.622149Figure
pd<- position_dodge(0.2)
g21338_fig<-ggplot(data=g21338_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(CUB~and~peptidase~domain-containing~protein~2-like~expression))+
ggtitle(~cyan)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg21338_fig
##L-type calcium channel alpha-1 subunit
wgcna_counts_filtered_long_g21501 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g21501.t1")
wgcna_counts_filtered_long_g21501 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 28 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 29 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 19 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 16 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 31 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 28 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 0 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 23 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 0 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 0 18 D Stable
## # ℹ 36 more rowsg21501.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g21501 , na.action=na.exclude)
car::Anova(g21501.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 33.2008 1 8.312e-09 ***
## Origin 0.5237 1 0.4693
## Treatment2 0.4276 1 0.5132
## Origin:Treatment2 0.1305 1 0.7179
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g21501.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 20.3 2.71 19 14.6 26.0
## RS 17.8 2.60 19 12.4 23.3
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 2.48 3.75 23 0.662 0.5145
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g21501_sum<-summarySE(wgcna_counts_filtered_long_g21501 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g21501_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 22.09091 10.02452 3.022505 6.734561
## 2 RF Variable 11 18.54545 12.91792 3.894900 8.678379
## 3 RS Stable 12 18.25000 14.12364 4.077144 8.973734
## 4 RS Variable 12 17.41667 13.22160 3.816746 8.400601Figure
pd<- position_dodge(0.2)
g21501_fig<-ggplot(data=g21501_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(L-type~calcium~channel~alpha-1~subunit~expression))+
ggtitle(~turquoise)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg21501_fig
##.Protocadherin
wgcna_counts_filtered_long_g21501 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g22388.t1")
wgcna_counts_filtered_long_g21501 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 133 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 138 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 163 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 156 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 164 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 181 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 130 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 185 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 157 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 144 18 D Stable
## # ℹ 36 more rowsg21501.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g21501 , na.action=na.exclude)
car::Anova(g21501.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 302.6581 1 <2e-16 ***
## Origin 1.4899 1 0.2222
## Treatment2 0.7306 1 0.3927
## Origin:Treatment2 0.3360 1 0.5621
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g21501.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 162 6.38 19 149 176
## RS 142 6.11 19 129 155
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 20.4 8.83 23 2.306 0.0305
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g21501_sum<-summarySE(wgcna_counts_filtered_long_g21501 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g21501_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 157.0000 20.72679 6.249364 13.92445
## 2 RF Variable 11 167.9091 28.94980 8.728693 19.44874
## 3 RS Stable 12 141.7500 37.05064 10.695599 23.54085
## 4 RS Variable 12 142.4167 29.92250 8.637882 19.01185Figure
pd<- position_dodge(0.2)
g21501_fig<-ggplot(data=g21501_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g21501 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Protocadherin~expression), limits=c(100,200))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg21501_fig## Warning: Removed 5 rows containing missing values (`geom_point()`).
##Acropora yongei Na+/Ca2+ exchanger
wgcna_counts_filtered_long_g24639 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g24639.t1")
wgcna_counts_filtered_long_g24639 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 70 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 79 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 66 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 108 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 58 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 94 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 90 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 60 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 62 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 62 18 D Stable
## # ℹ 36 more rowsg24639.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g24639 , na.action=na.exclude)
car::Anova(g24639.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 208.9868 1 <2e-16 ***
## Origin 0.0585 1 0.8089
## Treatment2 0.0158 1 0.8999
## Origin:Treatment2 0.6123 1 0.4339
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g24639.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 80.6 4.21 19 71.8 89.4
## RS 78.5 4.03 19 70.1 86.9
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 2.1 5.6 23 0.375 0.7115
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g24639_sum<-summarySE(wgcna_counts_filtered_long_g24639 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g24639_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 81.27273 15.07376 4.544909 10.12669
## 2 RF Variable 11 82.18182 20.54175 6.193572 13.80014
## 3 RS Stable 12 81.58333 16.05365 4.634290 10.20000
## 4 RS Variable 12 74.66667 21.21035 6.122900 13.47641Figure
pd<- position_dodge(0.2)
g24639_fig<-ggplot(data=g24639_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Acropora~yongei~"Na+/Ca2+"~exchanger~expression))+
ggtitle(~tan)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg24639_fig
##**mammalian ependymin-related protein 1-like
wgcna_counts_filtered_long_g25351 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g25351.t1")
wgcna_counts_filtered_long_g25351 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 5677 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 5022 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 13830 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 5729 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 5012 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 8437 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 7397 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 6199 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 5018 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 4565 18 D Stable
## # ℹ 36 more rowsg25351.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g25351 , na.action=na.exclude)
car::Anova(g25351.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 113.2918 1 < 2e-16 ***
## Origin 6.5397 1 0.01055 *
## Treatment2 2.9203 1 0.08747 .
## Origin:Treatment2 1.6792 1 0.19503
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g25351.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 6638 396 19 5809 7466
## RS 3944 379 19 3151 4738
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 2693 548 23 4.912 0.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g25351_sum<-summarySE(wgcna_counts_filtered_long_g25351 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g25351_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 5960.909 1716.8516 517.6502 1153.3966
## 2 RF Variable 11 7314.364 2894.5978 872.7541 1944.6172
## 3 RS Stable 12 3978.167 792.8993 228.8903 503.7842
## 4 RS Variable 12 3910.750 1499.1418 432.7650 952.5093Figure
pd<- position_dodge(0.2)
g25351_fig<-ggplot(data=g25351_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g25351,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(mammalian~ependymin-related~protein~1-like~expression), limits=c(2500,12500))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg25351_fig## Warning: Removed 2 rows containing missing values (`geom_point()`).
##MAM and LDLr domain-containing protein
wgcna_counts_filtered_long_g25935 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g25935.t1")
wgcna_counts_filtered_long_g25935 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 1097 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 929 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 1262 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 1089 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 1127 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 1375 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 1184 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 1275 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 1207 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 1525 18 D Stable
## # ℹ 36 more rowsg25935.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g25935 , na.action=na.exclude)
car::Anova(g25935.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 274.9500 1 <2e-16 ***
## Origin 2.0555 1 0.1517
## Treatment2 0.1990 1 0.6555
## Origin:Treatment2 0.4517 1 0.5015
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g25935.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 1324 59.6 19 1199 1449
## RS 1216 57.1 19 1096 1335
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 108 81 23 1.338 0.1941
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g25935_sum<-summarySE(wgcna_counts_filtered_long_g25935 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g25935_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 1344.182 302.1115 91.09005 202.9613
## 2 RF Variable 11 1294.818 214.8175 64.76992 144.3164
## 3 RS Stable 12 1183.250 328.6159 94.86323 208.7926
## 4 RS Variable 12 1236.833 214.0076 61.77867 135.9739Figure
pd<- position_dodge(0.2)
g25935_fig<-ggplot(data=g25935_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(MAM~and~LDLr~domain-containing~protein~expression))+
ggtitle(~brown)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg25935_fig
#SLIT-ROBO Rho GTPase-activating protein 1-like
wgcna_counts_filtered_long_g27376 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g27376.t1")
wgcna_counts_filtered_long_g27376 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 35 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 28 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 40 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 67 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 19 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 38 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 31 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 39 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 38 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 26 18 D Stable
## # ℹ 36 more rowsg27376.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g27376 , na.action=na.exclude)
car::Anova(g27376.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 98.9685 1 <2e-16 ***
## Origin 0.0141 1 0.9054
## Treatment2 0.2753 1 0.5998
## Origin:Treatment2 0.7895 1 0.3742
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g27376.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 36.0 2.72 19 30.4 41.7
## RS 38.6 2.60 19 33.2 44.1
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -2.58 3.72 23 -0.694 0.4945
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g27376_sum<-summarySE(wgcna_counts_filtered_long_g27376 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g27376_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 37.36364 10.763575 3.245340 7.231068
## 2 RF Variable 11 34.63636 11.209574 3.379814 7.530694
## 3 RS Stable 12 36.83333 16.889121 4.875469 10.730836
## 4 RS Variable 12 40.50000 9.491623 2.739996 6.030690Figure
pd<- position_dodge(0.2)
g27376_fig<-ggplot(data=g27376_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(SLIT-ROBO~Rho~GTPase-activating~protein~1-like~expression))+
ggtitle(~grey60)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg27376_fig
##**Hephaestin-like protein - frontloaded
wgcna_counts_filtered_long_g27566 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g27566.t1")
wgcna_counts_filtered_long_g27566 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 12 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 24 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 15 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 21 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 30 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 26 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 15 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 15 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 0 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 17 18 D Stable
## # ℹ 36 more rowsg27566.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g27566 , na.action=na.exclude)
car::Anova(g27566.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 56.1260 1 6.797e-14 ***
## Origin 6.3687 1 0.01162 *
## Treatment2 0.0303 1 0.86179
## Origin:Treatment2 0.0252 1 0.87392
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g27566.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 19.0 1.83 19 15.22 22.9
## RS 10.4 1.75 19 6.75 14.1
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 8.63 2.53 23 3.410 0.0024
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g27566_sum<-summarySE(wgcna_counts_filtered_long_g27566 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g27566_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 19.36364 8.139690 2.454209 5.468318
## 2 RF Variable 11 18.72727 8.978763 2.707199 6.032015
## 3 RS Stable 12 10.33333 9.528267 2.750574 6.053972
## 4 RS Variable 12 10.50000 7.501515 2.165501 4.766235Figure
pd<- position_dodge(0.2)
g27566_fig<-ggplot(data=g27566_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g27566 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Hephaestin-like~protein~expression))+
ggtitle(~brown)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg27566_fig
plasma membrane calcium ATPase - frontloaded
wgcna_counts_filtered_long_g27976 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g27976.t1")
wgcna_counts_filtered_long_g27976 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 197 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 219 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 204 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 124 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 174 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 155 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 183 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 173 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 122 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 223 18 D Stable
## # ℹ 36 more rowsg27976.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g27976 , na.action=na.exclude)
car::Anova(g27976.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 209.0712 1 <2e-16 ***
## Origin 0.0008 1 0.9774
## Treatment2 0.3269 1 0.5675
## Origin:Treatment2 0.0012 1 0.9725
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g27976.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 184 8.74 19 166 202
## RS 184 8.37 19 166 201
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 0.0682 12.1 23 0.006 0.9956
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g27976_sum<-summarySE(wgcna_counts_filtered_long_g27976 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g27976_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 178.8182 39.85177 12.01576 26.77278
## 2 RF Variable 11 188.8182 33.30111 10.04066 22.37199
## 3 RS Stable 12 178.3333 45.44794 13.11969 28.87624
## 4 RS Variable 12 189.1667 43.65950 12.60341 27.73992Figure
pd<- position_dodge(0.2)
g27976_fig<-ggplot(data=g27976_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(plasma~membrane~calcium~ATPase~expression))+
ggtitle(~brown)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg27976_fig
#von Willebrand factor D and EGF domain-containing protein-like
wgcna_counts_filtered_long_g28226 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g28226.t2")
wgcna_counts_filtered_long_g28226 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 0 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 0 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 0 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 0 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 66 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 103 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 0 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 0 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 0 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 0 18 D Stable
## # ℹ 36 more rowsg28226.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g28226 , na.action=na.exclude)
car::Anova(g28226.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 13.4863 1 0.0002403 ***
## Origin 6.2703 1 0.0122776 *
## Treatment2 2.3220 1 0.1275576
## Origin:Treatment2 3.0108 1 0.0827128 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g28226.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 30.2 8.97 19 11.43 49.0
## RS 16.6 8.82 19 -1.88 35.1
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 13.6 7.17 23 1.899 0.0702
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g28226_sum<-summarySE(wgcna_counts_filtered_long_g28226 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g28226_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 35.18182 57.46968 17.327759 38.60865
## 2 RF Variable 11 26.36364 40.63317 12.251362 27.29774
## 3 RS Stable 12 13.08333 16.41761 4.739355 10.43125
## 4 RS Variable 12 18.16667 25.95392 7.492252 16.49034Figure
pd<- position_dodge(0.2)
g28226_fig<-ggplot(data=g28226_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g28226,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(von~Willebrand~factor~D~and~EGF~domain-containing~protein-like~expression), limits=c(0,120))+
ggtitle(~magenta)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg28226_fig## Warning: Removed 1 rows containing missing values (`geom_point()`).
#low-density lipoprotein receptor-related protein 8-like
wgcna_counts_filtered_long_g4085 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g4085.t1")
wgcna_counts_filtered_long_g4085 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 37 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 42 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 54 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 47 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 54 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 73 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 40 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 31 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 52 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 49 18 D Stable
## # ℹ 36 more rowsg4085.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g4085 , na.action=na.exclude)
car::Anova(g4085.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 247.3507 1 <2e-16 ***
## Origin 0.9727 1 0.3240
## Treatment2 0.8654 1 0.3522
## Origin:Treatment2 2.5023 1 0.1137
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g4085.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 45.8 2.15 19 41.3 50.3
## RS 46.4 2.06 19 42.1 50.7
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -0.557 2.98 23 -0.187 0.8532
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g4085_sum<-summarySE(wgcna_counts_filtered_long_g4085 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g4085_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 47.81818 11.151845 3.362408 7.491911
## 2 RF Variable 11 43.81818 8.364427 2.521970 5.619298
## 3 RS Stable 12 43.66667 8.391915 2.422537 5.331969
## 4 RS Variable 12 49.08333 11.881677 3.429945 7.549257Figure
pd<- position_dodge(0.2)
g4085_fig<-ggplot(data=g4085_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(low-density~lipoprotein~receptor-related~protein~8-like~expression))+
ggtitle(~cyan)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg4085_fig
#**Cephalotoxin-like protein
wgcna_counts_filtered_long_g5013 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g5013.t1")
wgcna_counts_filtered_long_g5013 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 6 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 14 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 14 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 6 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 23 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 9 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 16 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 21 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 22 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 12 18 D Stable
## # ℹ 36 more rowsg5013.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g5013 , na.action=na.exclude)
car::Anova(g5013.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 44.8432 1 2.135e-11 ***
## Origin 6.0215 1 0.01413 *
## Treatment2 0.2947 1 0.58719
## Origin:Treatment2 0.0002 1 0.99022
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g5013.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 13.76 1.41 19 10.80 16.7
## RS 7.22 1.35 19 4.39 10.0
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 6.54 1.93 23 3.393 0.0025
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g5013_sum<-summarySE(wgcna_counts_filtered_long_g5013 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g5013_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 13.090909 7.063350 2.129680 4.745223
## 2 RF Variable 11 14.545455 6.170310 1.860419 4.145271
## 3 RS Stable 12 6.416667 5.583390 1.611786 3.547517
## 4 RS Variable 12 7.916667 6.934215 2.001735 4.405790Figure
pd<- position_dodge(0.2)
g5013_fig<-ggplot(data=g5013_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g5013 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(cephalotoxin-like~expression))+
ggtitle(~blue)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg5013_fig
#**sodium bicarbonate cotransporter 3-like isoform X - SLC4A7
wgcna_counts_filtered_long_g7402 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g7402.t1")
wgcna_counts_filtered_long_g7402 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 391 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 785 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 328 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 522 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 799 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 863 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 557 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 533 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 906 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 818 18 D Stable
## # ℹ 36 more rowsg7402.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g7402, na.action=na.exclude)
car::Anova(g7402.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 102.5059 1 < 2e-16 ***
## Origin 6.4414 1 0.01115 *
## Treatment2 1.0268 1 0.31091
## Origin:Treatment2 0.6016 1 0.43798
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g7402.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 660 49.6 19 556 764
## RS 467 47.5 19 367 566
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 193 68.7 23 2.814 0.0099
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g7402_sum<-summarySE(wgcna_counts_filtered_long_g7402 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g7402_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 710.3636 265.8632 80.16078 178.6093
## 2 RF Variable 11 609.8182 241.3930 72.78272 162.1700
## 3 RS Stable 12 463.8333 248.6454 71.77773 157.9817
## 4 RS Variable 12 469.8333 166.4407 48.04730 105.7514Figure
pd<- position_dodge(0.2)
g7402_fig<-ggplot(data=g7402_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g7402 ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(sodium~bicarbonate~cotransporter~3-like~isoform~X2~expression))+
ggtitle(~pink)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg7402_fig
##*protein lingerer-like
wgcna_counts_filtered_long_g7902 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___RNAseq.g7908.t1")
wgcna_counts_filtered_long_g7902 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 51 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 73 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 101 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 79 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 107 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 95 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 61 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 75 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 92 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 68 18 D Stable
## # ℹ 36 more rowsg7902.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g7902, na.action=na.exclude)
car::Anova(g7902.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 167.9194 1 < 2.2e-16 ***
## Origin 11.0673 1 0.0008786 ***
## Treatment2 0.0009 1 0.9759807
## Origin:Treatment2 0.0806 1 0.7764657
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g7902.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 82.9 4.53 19 73.4 92.3
## RS 51.6 4.34 19 42.5 60.7
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 31.3 6.27 23 4.989 <.0001
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g7902_sum<-summarySE(wgcna_counts_filtered_long_g7902 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g7902_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 83.00000 24.55606 7.403930 16.49698
## 2 RF Variable 11 82.72727 21.99587 6.632004 14.77702
## 3 RS Stable 12 53.50000 19.52853 5.637402 12.40784
## 4 RS Variable 12 49.66667 18.80683 5.429065 11.94929Figure
pd<- position_dodge(0.2)
g7902_fig<-ggplot(data=g7902_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_g7902,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(protein~lingerer-like~expression))+
ggtitle(~pink)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg7902_fig
##***CA1
wgcna_counts_filtered_long_CA1<- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___TS.g12304.t1")
wgcna_counts_filtered_long_CA1## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 2854 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 4635 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 2949 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 4681 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 7704 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 8665 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 3948 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 3887 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 6896 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 6597 18 D Stable
## # ℹ 36 more rowsCA1.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_CA1, na.action=na.exclude)
car::Anova(CA1.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 97.1566 1 < 2.2e-16 ***
## Origin 15.7334 1 7.293e-05 ***
## Treatment2 2.1565 1 0.142
## Origin:Treatment2 1.3573 1 0.244
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(CA1.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 5346 434 19 4437 6256
## RS 2716 416 19 1846 3587
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 2630 597 23 4.408 0.0002
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
CA1_sum<-summarySE(wgcna_counts_filtered_long_CA1, measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
CA1_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 5970.636 2621.084 790.2864 1760.8679
## 2 RF Variable 11 4733.455 1986.062 598.8204 1334.2549
## 3 RS Stable 12 2653.583 1847.177 533.2340 1173.6400
## 4 RS Variable 12 2775.250 1463.636 422.5154 929.9501Figure
pd<- position_dodge(0.2)
CA1_fig<-ggplot(data=CA1_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
geom_point(data=wgcna_counts_filtered_long_CA1,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(CA1~expression))+
ggtitle(~pink)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statCA1_fig
##Uncharacterized skeletal organic matrix protein-6 (USOMP6)
wgcna_counts_filtered_long_g22622 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___TS.g22622.t1")
wgcna_counts_filtered_long_g22622 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 902 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 733 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 1136 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 857 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 558 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 726 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 774 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 888 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 617 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 1117 18 D Stable
## # ℹ 36 more rowsg22622.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g22622, na.action=na.exclude)
car::Anova(g22622.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 155.9337 1 <2e-16 ***
## Origin 0.0013 1 0.9712
## Treatment2 0.5378 1 0.4633
## Origin:Treatment2 0.7739 1 0.3790
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g22622.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 903 54.2 19 790 1016
## RS 846 51.9 19 738 955
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 56.9 71 23 0.801 0.4313
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g22622_sum<-summarySE(wgcna_counts_filtered_long_g22622 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g22622_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 870.9091 172.5407 52.02298 115.9144
## 2 RF Variable 11 935.3636 251.5207 75.83634 168.9739
## 3 RS Stable 12 859.9167 294.4059 84.98765 187.0566
## 4 RS Variable 12 817.3333 190.8728 55.10023 121.2748Figure
pd<- position_dodge(0.2)
g22622_fig<-ggplot(data=g22622_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Uncharacterized~skeletal~organic~matrix~protein-6~(USOMP6)~expression))+
ggtitle(~green)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg22622_fig
##band 3 anion transport protein-like
wgcna_counts_filtered_long_g27873 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___TS.g27873.t1")
wgcna_counts_filtered_long_g27873 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 52 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 65 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 105 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 52 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 22 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 50 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 44 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 94 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 58 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 60 18 D Stable
## # ℹ 36 more rowsg27873.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g27873, na.action=na.exclude)
car::Anova(g27873.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 109.0383 1 <2e-16 ***
## Origin 0.7687 1 0.3806
## Treatment2 0.0156 1 0.9006
## Origin:Treatment2 0.2677 1 0.6049
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g27873.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 71.1 5.15 19 60.3 81.8
## RS 59.9 4.94 19 49.6 70.3
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 11.1 6.84 23 1.628 0.1171
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g27873_sum<-summarySE(wgcna_counts_filtered_long_g27873 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g27873_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 70.63636 21.65767 6.530032 14.54982
## 2 RF Variable 11 71.72727 27.65896 8.339491 18.58154
## 3 RS Stable 12 61.50000 22.71763 6.558016 14.43410
## 4 RS Variable 12 56.33333 17.55166 5.066726 11.15179Figure
pd<- position_dodge(0.2)
g27873_fig<-ggplot(data=g27873_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(band~3~anion~transport~protein-like~expression))+
ggtitle(~green)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg27873_fig
##collagenase 3-like
wgcna_counts_filtered_long_g5338 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___TS.g5338.t1")
wgcna_counts_filtered_long_g5338 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 15 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 14 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 9 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 35 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 22 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 16 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 24 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 30 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 21 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 22 18 D Stable
## # ℹ 36 more rowsg5338.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g5338, na.action=na.exclude)
car::Anova(g5338.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 64.2427 1 1.1e-15 ***
## Origin 0.7901 1 0.3741
## Treatment2 0.0652 1 0.7985
## Origin:Treatment2 0.0281 1 0.8670
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g5338.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 20.6 1.78 19 16.9 24.4
## RS 17.1 1.70 19 13.6 20.7
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS 3.51 2.46 23 1.425 0.1677
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g5338_sum<-summarySE(wgcna_counts_filtered_long_g5338 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g5338_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 20.18182 7.180782 2.165087 4.824115
## 2 RF Variable 11 21.09091 7.955558 2.398691 5.344617
## 3 RS Stable 12 17.08333 10.202124 2.945100 6.482120
## 4 RS Variable 12 17.16667 7.601834 2.194460 4.829975Figure
pd<- position_dodge(0.2)
g5338_fig<-ggplot(data=g5338_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(collagenase~3-like~expression))+
ggtitle(~green)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg5338_fig
Protocadherin
wgcna_counts_filtered_long_g6583 <- wgcna_counts_filtered_long %>%
filter(Gene == "Pocillopora_acuta_HIv2___TS.g6583.t1")
wgcna_counts_filtered_long_g6583 ## # A tibble: 46 × 7
## Gene Origin Colony.number Counts Colony Treatment Treatment2
## <chr> <fct> <chr> <int> <dbl> <fct> <fct>
## 1 Pocillopora_acuta_HI… RF 13B 254 13 B Variable
## 2 Pocillopora_acuta_HI… RF 13D 287 13 D Stable
## 3 Pocillopora_acuta_HI… RF 14B 305 14 B Variable
## 4 Pocillopora_acuta_HI… RF 14C 229 14 C Stable
## 5 Pocillopora_acuta_HI… RF 15B 320 15 B Variable
## 6 Pocillopora_acuta_HI… RF 15D 241 15 D Stable
## 7 Pocillopora_acuta_HI… RF 17B 171 17 B Variable
## 8 Pocillopora_acuta_HI… RF 17D 261 17 D Stable
## 9 Pocillopora_acuta_HI… RF 18B 256 18 B Variable
## 10 Pocillopora_acuta_HI… RF 18D 274 18 D Stable
## # ℹ 36 more rowsg6583.lme <- lme(Counts~Origin*Treatment2, random = ~1|Colony, data=wgcna_counts_filtered_long_g6583, na.action=na.exclude)
car::Anova(g6583.lme, type=3)## Analysis of Deviance Table (Type III tests)
##
## Response: Counts
## Chisq Df Pr(>Chisq)
## (Intercept) 245.3182 1 <2e-16 ***
## Origin 1.5732 1 0.2097
## Treatment2 0.0685 1 0.7935
## Origin:Treatment2 0.5712 1 0.4498
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1tukey3<- emmeans(g6583.lme, list(pairwise ~ Origin), adjust = "tukey")## NOTE: Results may be misleading due to involvement in interactionstukey3## $`emmeans of Origin`
## Origin emmean SE df lower.CL upper.CL
## RF 283 13.2 19 256 311
## RS 302 12.6 19 275 328
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
##
## $`pairwise differences of Origin`
## 1 estimate SE df t.ratio p.value
## RF - RS -18.3 18.1 23 -1.014 0.3210
##
## Results are averaged over the levels of: Treatment2
## Degrees-of-freedom method: containmentlibrary(Rmisc)
g6583_sum<-summarySE(wgcna_counts_filtered_long_g6583 , measurevar='Counts', groupvars=c('Origin', 'Treatment2'), na.rm=TRUE, conf.interval = 0.95)
g6583_sum## Origin Treatment2 N Counts sd se ci
## 1 RF Stable 11 286.9091 43.36463 13.07493 29.13275
## 2 RF Variable 11 280.2727 48.05224 14.48830 32.28194
## 3 RS Stable 12 318.0000 80.29491 23.17914 51.01695
## 4 RS Variable 12 284.8333 61.77648 17.83333 39.25090Figure
pd<- position_dodge(0.2)
g6583_fig<-ggplot(data=g6583_sum, aes(y=Counts, x=Treatment2, color=Origin),group = interaction(Origin))+
#geom_point(data=wgcna_counts_filtered_long_PRKCD ,aes(y=Counts, x=Treatment2, color=Origin), alpha=0.4, position = pd)+
geom_line(aes(group = interaction(Origin), stat="identity"),position=position_dodge(0.2))+
geom_point(size=3, stat="identity", position = pd)+
geom_errorbar(aes(ymin=Counts-se, ymax=Counts+se), stat="identity",width=0.2, position = pd)+
scale_fill_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_color_manual("Origin", values=c("RS"='paleturquoise3', "RF"= "indianred"))+
scale_y_continuous(expression(Protocadherin~expression))+
ggtitle(~black)+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5,size=12),#angling the labels on the x-axis
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),#telling it where to position our plot title
panel.background= element_rect(fill=NA, color='black'),#this is making the black box around the graph
#strip.background = element_blank(),
#strip.text = element_blank(),
legend.title = element_text(vjust=0.5,size=12),
legend.position="none",
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_text(size=12))#making the axis title larger ## Warning in geom_line(aes(group = interaction(Origin), stat = "identity"), :
## Ignoring unknown aesthetics: statg6583_fig
#Comparison of sig genes
biomin_compare_figs<-cowplot::plot_grid(g10093_fig,CA2_fig,g25351_fig,g5013_fig,g15280_fig, CA1_fig,g7402_fig,g7902_fig,g27566_fig,g28226_fig,g14505_fig,g11609_fig,g10093_fig, nrow=4)## Warning: Removed 2 rows containing missing values (`geom_point()`).## Warning: Removed 1 rows containing missing values (`geom_point()`).biomin_compare_figs
Merge the frequency of GOterms for both up- and down-regulation of calcification
cal_freq_terms <- merge(result,result_down, by=c("ParentTerm","Number.of.terms","Calcification.direction"),all=T)cal_biomin_terms <-read.csv("../../output/Biomineralization_goterms.csv")
head(cal_biomin_terms)## X.1 X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 1 GO:0006325 0.009864096 1 2
## 2 2 2 GO:0016570 0.009864096 1 2
## 3 3 3 GO:0051276 0.009864096 1 2
## 4 4 4 GO:0060537 0.012664490 1 2
## 5 5 5 GO:0000278 0.012973357 1 2
## 6 6 6 GO:0007049 0.012973357 1 2
## numInCat term ontology bh_adjust
## 1 2 chromatin organization BP 0.8862339
## 2 2 histone modification BP 0.8862339
## 3 2 chromosome organization BP 0.8862339
## 4 2 muscle tissue development BP 0.8862339
## 5 2 mitotic cell cycle BP 0.8862339
## 6 2 cell cycle BP 0.8862339
## ParentTerm Factor
## 1 chromatin remodeling Biomin
## 2 macromolecule deacylation Biomin
## 3 chromatin remodeling Biomin
## 4 muscle tissue development Biomin
## 5 microtubule cytoskeleton organization involved in mitosis Biomin
## 6 microtubule cytoskeleton organization involved in mitosis Biomincal_up_terms <-read.csv("../../output/WGCNA/GO_analysis/goseq_pattern_calcification_filtered.csv")
cal_up_terms<- cal_up_terms %>%
mutate(Factor = "Up")
head(cal_up_terms)## X.1 X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 1 GO:0000003 0 1 465
## 2 2 2 GO:0006139 0 1 608
## 3 3 3 GO:0006355 0 1 462
## 4 4 4 GO:0006725 0 1 650
## 5 5 5 GO:0006807 0 1 1215
## 6 6 6 GO:0006810 0 1 664
## numInCat term ontology bh_adjust
## 1 465 reproduction BP 0
## 2 608 nucleobase-containing compound metabolic process BP 0
## 3 462 regulation of DNA-templated transcription BP 0
## 4 650 cellular aromatic compound metabolic process BP 0
## 5 1215 nitrogen compound metabolic process BP 0
## 6 664 transport BP 0
## ParentTerm Factor
## 1 reproduction Up
## 2 metabolic process Up
## 3 regulation of biosynthetic process Up
## 4 metabolic process Up
## 5 metabolic process Up
## 6 localization Upcal_down_terms <-read.csv("../../output/WGCNA/GO_analysis/goseq_pattern_calcification_down_filtered.csv")
cal_down_terms<-cal_down_terms %>%
mutate(Factor = "Down")
head(cal_down_terms)## X.1 X GOterm over_represented_pvalue under_represented_pvalue numDEInCat
## 1 1 1 GO:0006139 0 1 485
## 2 2 2 GO:0006725 0 1 513
## 3 3 3 GO:0006807 0 1 910
## 4 4 4 GO:0006810 0 1 419
## 5 5 5 GO:0006950 0 1 364
## 6 6 6 GO:0006996 0 1 461
## numInCat term ontology bh_adjust
## 1 485 nucleobase-containing compound metabolic process BP 0
## 2 513 cellular aromatic compound metabolic process BP 0
## 3 910 nitrogen compound metabolic process BP 0
## 4 419 transport BP 0
## 5 364 response to stress BP 0
## 6 461 organelle organization BP 0
## ParentTerm Factor
## 1 gene expression Down
## 2 gene expression Down
## 3 gene expression Down
## 4 localization Down
## 5 response to stress Down
## 6 cellular component organization or biogenesis Downcolnames(cal_biomin_terms)## [1] "X.1" "X"
## [3] "GOterm" "over_represented_pvalue"
## [5] "under_represented_pvalue" "numDEInCat"
## [7] "numInCat" "term"
## [9] "ontology" "bh_adjust"
## [11] "ParentTerm" "Factor"colnames(cal_up_terms)## [1] "X.1" "X"
## [3] "GOterm" "over_represented_pvalue"
## [5] "under_represented_pvalue" "numDEInCat"
## [7] "numInCat" "term"
## [9] "ontology" "bh_adjust"
## [11] "ParentTerm" "Factor"colnames(cal_down_terms)## [1] "X.1" "X"
## [3] "GOterm" "over_represented_pvalue"
## [5] "under_represented_pvalue" "numDEInCat"
## [7] "numInCat" "term"
## [9] "ontology" "bh_adjust"
## [11] "ParentTerm" "Factor"colnames(cal_freq_terms)## [1] "ParentTerm" "Number.of.terms"
## [3] "Calcification.direction"Merge biomineralization, up and down-regulation of calcification GOterms
all_terms<- merge(cal_up_terms,cal_down_terms, by=c("Factor","GOterm","X.1","X","GOterm","over_represented_pvalue","under_represented_pvalue","numDEInCat","numInCat","term","ontology","bh_adjust","ParentTerm"),all=T)
all_terms<- merge(all_terms,cal_biomin_terms, by=c("Factor","GOterm","X.1","X","GOterm","over_represented_pvalue","under_represented_pvalue","numDEInCat","numInCat","term","ontology","bh_adjust","ParentTerm"),all=T)
all_terms$GOterm<-as.factor(all_terms$GOterm)
head(all_terms)## Factor GOterm X.1 X over_represented_pvalue under_represented_pvalue
## 1 Biomin GO:0000003 300 343 0.54128075 0.8456150
## 2 Biomin GO:0000041 417 565 1.00000000 0.8535542
## 3 Biomin GO:0000122 65 70 0.12290664 1.0000000
## 4 Biomin GO:0000132 107 122 0.15108190 1.0000000
## 5 Biomin GO:0000226 256 298 0.37841205 0.9418666
## 6 Biomin GO:0000278 5 5 0.01297336 1.0000000
## numDEInCat numInCat term
## 1 1 5 reproduction
## 2 0 1 transition metal ion transport
## 3 1 1 negative regulation of transcription by RNA polymerase II
## 4 1 1 establishment of mitotic spindle orientation
## 5 1 3 microtubule cytoskeleton organization
## 6 2 2 mitotic cell cycle
## ontology bh_adjust ParentTerm
## 1 BP 1.0000000 female sex differentiation
## 2 BP 1.0000000 calcium ion import
## 3 BP 0.8862339 negative regulation of biological process
## 4 BP 0.8862339 establishment of mitotic spindle orientation
## 5 BP 1.0000000 microtubule-based process
## 6 BP 0.8862339 microtubule cytoskeleton organization involved in mitosistail(all_terms)## Factor GOterm X.1 X over_represented_pvalue
## 4614 Up GO:2001242 922 983 5.023572e-31
## 4615 Up GO:2001243 1342 1500 2.171774e-18
## 4616 Up GO:2001251 978 1045 3.318296e-29
## 4617 Up GO:2001252 713 759 1.527374e-42
## 4618 Up GO:2001257 1514 1737 2.346995e-15
## 4619 Up GO:2001259 2050 2535 2.136874e-09
## under_represented_pvalue numDEInCat numInCat
## 4614 1 42 42
## 4615 1 24 24
## 4616 1 43 43
## 4617 1 61 61
## 4618 1 21 21
## 4619 1 12 12
## term ontology
## 4614 regulation of intrinsic apoptotic signaling pathway BP
## 4615 negative regulation of intrinsic apoptotic signaling pathway BP
## 4616 negative regulation of chromosome organization BP
## 4617 positive regulation of chromosome organization BP
## 4618 regulation of cation channel activity BP
## 4619 positive regulation of cation channel activity BP
## bh_adjust ParentTerm
## 4614 5.462041e-30 regulation of cell death
## 4615 1.558880e-17 regulation of cell death
## 4616 3.389807e-28 regulation of cellular component organization
## 4617 2.152863e-41 regulation of cellular component organization
## 4618 1.467562e-14 regulation of localization
## 4619 9.175464e-09 regulation of localizationgoterms_shared <- all_terms %>%
group_by(GOterm) %>%
dplyr::summarise(
ParentTerm = paste(unique(ParentTerm), collapse = ", "),
Factor = paste(unique(Factor), collapse = ", ")
)
dim(goterms_shared)## [1] 2322 3write.csv(goterms_shared, "../../output/WGCNA/GO_analysis/Merged_GOterms_factor_ParentTerm.csv")result_unique <- goterms_shared %>%
group_by(ParentTerm,Factor) %>%
dplyr::summarise(SharedGOterms = n_distinct(GOterm))%>%
arrange(-SharedGOterms)## `summarise()` has grouped output by 'ParentTerm'. You can override using the
## `.groups` argument.result_filtered_up<- result_unique %>%
dplyr::filter(Factor=="Biomin, Down, Up" | Factor=="Biomin, Up")
#dplyr::filter(SharedGOterms>=5) result_filtered_down<- result_unique %>%
dplyr::filter(Factor=="Biomin, Down, Up" | Factor=="Biomin, Down")
#dplyr::filter(SharedGOterms>=5) Calculate the percentage of shared GOterms downregulation of calcification
merged_down <- result_filtered_down %>%
tidyr::separate_rows(ParentTerm, sep = ", ") %>%
dplyr::group_by(ParentTerm) %>%
dplyr::summarize(Sum_of_SharedGOterms = sum(SharedGOterms, na.rm = TRUE))%>%
left_join(result_down, by = "ParentTerm") %>%
mutate(Proportion.of.shared.GO.terms.with.biomineralization.genes = Sum_of_SharedGOterms / Number.of.terms)%>%
mutate(Percentage.of.shared.GO.terms.with.biomineralization.genes = Proportion.of.shared.GO.terms.with.biomineralization.genes * 100)
merged_down_clean<- na.omit(merged_down)
merged_down_clean## # A tibble: 96 × 6
## ParentTerm Sum_of_SharedGOterms Number.of.terms Calcification.direct…¹
## <chr> <int> <int> <chr>
## 1 actin filament-b… 5 11 Down
## 2 aging 1 2 Down
## 3 amide metabolic … 2 12 Down
## 4 ammonium ion met… 1 1 Down
## 5 anatomical struc… 23 29 Down
## 6 animal organ dev… 12 16 Down
## 7 behavior 3 8 Down
## 8 biosynthetic pro… 1 24 Down
## 9 carbohydrate der… 4 18 Down
## 10 catabolic process 15 38 Down
## # ℹ 86 more rows
## # ℹ abbreviated name: ¹​Calcification.direction
## # ℹ 2 more variables:
## # Proportion.of.shared.GO.terms.with.biomineralization.genes <dbl>,
## # Percentage.of.shared.GO.terms.with.biomineralization.genes <dbl>###Frequency >10 up
cal_freq_terms_filtered_up <- merged_up_clean %>%
filter(Number.of.terms>=10) %>%
filter(Calcification.direction=="Up")
cal_freq_terms_filtered_up## # A tibble: 71 × 6
## ParentTerm Sum_of_SharedGOterms Number.of.terms Calcification.direct…¹
## <chr> <int> <int> <chr>
## 1 actin filament-b… 5 12 Up
## 2 amide metabolic … 2 12 Up
## 3 anatomical struc… 24 27 Up
## 4 animal organ dev… 12 15 Up
## 5 behavior 3 10 Up
## 6 biosynthetic pro… 1 25 Up
## 7 carbohydrate der… 4 27 Up
## 8 catabolic process 15 39 Up
## 9 cell cycle 9 63 Up
## 10 cell division 4 14 Up
## # ℹ 61 more rows
## # ℹ abbreviated name: ¹​Calcification.direction
## # ℹ 2 more variables:
## # Proportion.of.shared.GO.terms.with.biomineralization.genes <dbl>,
## # Percentage.of.shared.GO.terms.with.biomineralization.genes <dbl>Figure
#counts$Direction.of.flat.origin<- factor(counts$Direction.of.flat.origin, levels =c("up","no pattern","down"))
#counts$Module<- factor(counts$Module, levels=c("Blue","Brown","Greenyellow","Cyan","Pink","Magenta","Lightcyan","Midnight blue","Purple","Turquiose","Red","Black"))
freq_fig_up<-ggplot(cal_freq_terms_filtered_up, aes(y=Number.of.terms,x=reorder(ParentTerm, Number.of.terms), fill=Percentage.of.shared.GO.terms.with.biomineralization.genes,group=1))+
#facet_wrap(~Calcification.direction, nrow = 1)+
geom_point(size=5, alpha=1, pch=21,color="black")+
geom_segment(aes(x=ParentTerm, xend=ParentTerm, y=0, yend=Number.of.terms)) +
geom_hline(yintercept = 0, linetype="solid", color = 'black', size=0.5, show.legend = TRUE)+
coord_flip()+
scale_y_continuous(expression(GO~term~counts),limits=c(0,70))+
#scale_color_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
#scale_fill_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
scale_fill_gradientn(colours=c("white","#fddbc7","#f4a582","#d6604d","#b2182b"), na.value = "grey98",limits = c(0, 100))+
#scale_color_gradientn(colours=c("#b2182b","#fddbc7","white","#d1e5f0","#67a9cf", "#67a9cf", "#2166ac"), na.value = "grey98",limits = c(-0, 40))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, hjust=0.95,size=12),
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),
panel.background= element_rect(fill=NA, color='black'),
legend.title = element_blank(),
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_blank(),
strip.text = element_text(size=12))#making the axis title larger ## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.freq_fig_up
ggsave(filename="../../output/WGCNA/GO_analysis/freq_fig_up.png", plot=freq_fig_up, dpi=300, height=12, units="in", limitsize=FALSE)## Saving 7 x 12 in image###Frequency >10 down
cal_freq_terms_filtered_down <- merged_down_clean %>%
filter(Number.of.terms>=10)
#filter(Calcification.direction=="Down")
cal_freq_terms_filtered_down## # A tibble: 56 × 6
## ParentTerm Sum_of_SharedGOterms Number.of.terms Calcification.direct…¹
## <chr> <int> <int> <chr>
## 1 actin filament-b… 5 11 Down
## 2 amide metabolic … 2 12 Down
## 3 anatomical struc… 23 29 Down
## 4 animal organ dev… 12 16 Down
## 5 biosynthetic pro… 1 24 Down
## 6 carbohydrate der… 4 18 Down
## 7 catabolic process 15 38 Down
## 8 cell cycle 9 43 Down
## 9 cell differentia… 11 12 Down
## 10 cell projection … 15 20 Down
## # ℹ 46 more rows
## # ℹ abbreviated name: ¹​Calcification.direction
## # ℹ 2 more variables:
## # Proportion.of.shared.GO.terms.with.biomineralization.genes <dbl>,
## # Percentage.of.shared.GO.terms.with.biomineralization.genes <dbl>Figure
freq_fig_down<-ggplot(cal_freq_terms_filtered_down, aes(y=Number.of.terms,x=reorder(ParentTerm, Number.of.terms), fill=Percentage.of.shared.GO.terms.with.biomineralization.genes))+
#facet_wrap(~Calcification.direction, nrow = 1)+
geom_point(size=5, alpha=1, pch=21,color="black")+
geom_segment(aes(x=ParentTerm, xend=ParentTerm, y=0, yend=Number.of.terms)) +
#geom_hline(yintercept = 0, linetype="solid", color = 'black', size=0.5, show.legend = TRUE)+
coord_flip()+
scale_y_continuous(expression(GO~term~counts),limits=c(0,70))+
#scale_color_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
#scale_fill_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
scale_fill_gradientn(colours=c("white","#d1e5f0","#92c5de","#4393c3","#2166ac"), na.value = "grey98",limits = c(0, 100))+
#scale_fill_gradientn(colours=c("#fddbc7","#f4a582","#d6604d","#b2182b"), na.value = "grey98",limits = c(10, 40))+
#scale_color_gradientn(colours=c("#b2182b","#fddbc7","white","#d1e5f0","#67a9cf", "#67a9cf", "#2166ac"), na.value = "grey98",limits = c(-0, 40))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, hjust=0.95, size=12),
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),
panel.background= element_rect(fill=NA, color='black'),
legend.title = element_blank(),
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_blank(),
strip.text = element_text(size=12))#making the axis title larger
freq_fig_down
ggsave(filename="../../output/WGCNA/GO_analysis/freq_fig_down.png", plot=freq_fig_down, dpi=300, height=12, units="in", limitsize=FALSE)## Saving 7 x 12 in imagecompare_figs<-cowplot::plot_grid(freq_fig_up, freq_fig_down, nrow=2, align="v")
compare_figs
All terms up
cal_freq_terms_filtered_up_all <- merged_up_clean %>%
filter(Calcification.direction=="Up")
cal_freq_terms_filtered_up_all## # A tibble: 106 × 6
## ParentTerm Sum_of_SharedGOterms Number.of.terms Calcification.direct…¹
## <chr> <int> <int> <chr>
## 1 actin filament-b… 5 12 Up
## 2 aging 1 2 Up
## 3 amide metabolic … 2 12 Up
## 4 ammonium ion met… 1 1 Up
## 5 anatomical struc… 24 27 Up
## 6 animal organ dev… 12 15 Up
## 7 behavior 3 10 Up
## 8 biosynthetic pro… 1 25 Up
## 9 carbohydrate der… 4 27 Up
## 10 catabolic process 15 39 Up
## # ℹ 96 more rows
## # ℹ abbreviated name: ¹​Calcification.direction
## # ℹ 2 more variables:
## # Proportion.of.shared.GO.terms.with.biomineralization.genes <dbl>,
## # Percentage.of.shared.GO.terms.with.biomineralization.genes <dbl>Figure
#counts$Direction.of.flat.origin<- factor(counts$Direction.of.flat.origin, levels =c("up","no pattern","down"))
#counts$Module<- factor(counts$Module, levels=c("Blue","Brown","Greenyellow","Cyan","Pink","Magenta","Lightcyan","Midnight blue","Purple","Turquiose","Red","Black"))
freq_fig_up<-ggplot(cal_freq_terms_filtered_up_all, aes(y=Number.of.terms,x=reorder(ParentTerm, Number.of.terms), fill=Percentage.of.shared.GO.terms.with.biomineralization.genes,group=1))+
#facet_wrap(~Calcification.direction, nrow = 1)+
geom_point(size=5, alpha=1, pch=21,color="black")+
geom_segment(aes(x=ParentTerm, xend=ParentTerm, y=0, yend=Number.of.terms)) +
geom_hline(yintercept = 0, linetype="solid", color = 'black', size=0.5, show.legend = TRUE)+
coord_flip()+
scale_y_continuous(expression(GO~term~counts),limits=c(0,70))+
#scale_color_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
#scale_fill_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
scale_fill_gradientn(colours=c("white","#fddbc7","#f4a582","#d6604d","#b2182b"), na.value = "grey98",limits = c(0, 100))+
#scale_color_gradientn(colours=c("#b2182b","#fddbc7","white","#d1e5f0","#67a9cf", "#67a9cf", "#2166ac"), na.value = "grey98",limits = c(-0, 40))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, hjust=0.95,size=12),
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),
panel.background= element_rect(fill=NA, color='black'),
legend.title = element_blank(),
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_blank(),
strip.text = element_text(size=12))#making the axis title larger
freq_fig_up
ggsave(filename="../../output/WGCNA/GO_analysis/freq_fig_up_ALL.png", plot=freq_fig_up, dpi=300, height=17, units="in", limitsize=FALSE)## Saving 7 x 17 in imageAll terms down
cal_freq_terms_filtered_down_all <- merged_down_clean %>%
filter(Calcification.direction=="Down")
cal_freq_terms_filtered_down_all## # A tibble: 96 × 6
## ParentTerm Sum_of_SharedGOterms Number.of.terms Calcification.direct…¹
## <chr> <int> <int> <chr>
## 1 actin filament-b… 5 11 Down
## 2 aging 1 2 Down
## 3 amide metabolic … 2 12 Down
## 4 ammonium ion met… 1 1 Down
## 5 anatomical struc… 23 29 Down
## 6 animal organ dev… 12 16 Down
## 7 behavior 3 8 Down
## 8 biosynthetic pro… 1 24 Down
## 9 carbohydrate der… 4 18 Down
## 10 catabolic process 15 38 Down
## # ℹ 86 more rows
## # ℹ abbreviated name: ¹​Calcification.direction
## # ℹ 2 more variables:
## # Proportion.of.shared.GO.terms.with.biomineralization.genes <dbl>,
## # Percentage.of.shared.GO.terms.with.biomineralization.genes <dbl>Figure
freq_fig_down<-ggplot(cal_freq_terms_filtered_down_all, aes(y=Number.of.terms,x=reorder(ParentTerm, Number.of.terms), fill=Percentage.of.shared.GO.terms.with.biomineralization.genes))+
#facet_wrap(~Calcification.direction, nrow = 1)+
geom_point(size=5, alpha=1, pch=21,color="black")+
geom_segment(aes(x=ParentTerm, xend=ParentTerm, y=0, yend=Number.of.terms)) +
#geom_hline(yintercept = 0, linetype="solid", color = 'black', size=0.5, show.legend = TRUE)+
coord_flip()+
scale_y_continuous(expression(GO~term~counts),limits=c(0,70))+
#scale_color_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
#scale_fill_manual("Direction.of.flat.origin",values= c("up"="#b2182b","no pattern"="grey","down" ="#67a9cf"))+
scale_fill_gradientn(colours=c("white","#d1e5f0","#92c5de","#4393c3","#2166ac"), na.value = "grey98",limits = c(0, 100))+
#scale_fill_gradientn(colours=c("#fddbc7","#f4a582","#d6604d","#b2182b"), na.value = "grey98",limits = c(10, 40))+
#scale_color_gradientn(colours=c("#b2182b","#fddbc7","white","#d1e5f0","#67a9cf", "#67a9cf", "#2166ac"), na.value = "grey98",limits = c(-0, 40))+
theme_classic()+
theme(axis.text.x=element_text(vjust=0.5, hjust=0.95, size=12),
plot.title = element_text(margin = margin(t = 10, b = 10), hjust=0.5),
panel.background= element_rect(fill=NA, color='black'),
legend.title = element_blank(),
axis.text.y = element_text(vjust=0.5, size=12), #making the axis text larger
axis.title.x = element_blank(),#making the axis title larger
axis.title.y = element_blank(),
strip.text = element_text(size=12))#making the axis title larger
freq_fig_down
ggsave(filename="../../output/WGCNA/GO_analysis/freq_fig_down_ALL.png", plot=freq_fig_down, dpi=300, height=17, units="in", limitsize=FALSE)## Saving 7 x 17 in imagecompare_figs<-cowplot::plot_grid(freq_fig_up, freq_fig_down, ncol=2, align="h")
compare_figs